Raster reprojection OpenLayers 3

Question

I'm trying to reproject a WMS source in EPSG:3395 to EPSG:32661 in OpenLayers 3. I have proj4js and the projection definitions loaded but nothing shows up.

Curiously, not even the example at http://openlayers.org/en/v3.12.1/examples/reprojection-by-code.html seems to be working.

Answer 1

The error happens when a null value is processed by the proj4 code called by OpenLayers when an object containing x and y properties is expected. Depending on the version of proj4.js being used it can happen in any transform where one of the coordinates in invalid in the other projection but seems much likely less when tranforming to or from EPSG:4326 The OpenLayers examples try to prevent errors by setting extents on projections, but that is undesirable if you want to show a country in its national projection but also the show surrounding countries (albeit with some distortion), or include adjacent areas in UTM projections. I have used ol.proj.addCoordinateTransforms to ensure all transforms go via EPSG:4326 which eliminates most of the errors and added try/catch to override any remaining errors (although if using the latest version 2.5.0 of proj4.js the try/catch seems redundant, even that version produces errors transforming directly between other projections). However some conversions from EPSG:2163 to EPSG:4326 needed for ol.proj.getPointResolution() fail even with the latest proj4.js making it impossible to center or zoom out the map at some locations. Fortunately EPSG:4269 is predefined in proj4 and when used as an intermediate instead of EPSG:4326 it allows try/catch handling of errors involving EPSG:4326 (unless the reprojection is to or from EPSG:4269 when EPSG:4326 should be used!) Here's a rewrite of the OpenLayers example with all transforms made via EPSG:4269 and no restrictions on projection extents

function reprojectionErrorHandler (projections, opt_intermediate) {

var intermediate = opt_intermediate || 'EPSG:4269';

function transform (projA, projB) {

    return function (input, opt_output, opt_dimension) {
        var length = input.length;
        var dimension = opt_dimension !== undefined ? opt_dimension : 2;
        var output = opt_output !== undefined ? opt_output : new Array(length);
        var ll, point, i, j;
        try {
            for (i = 0; i < length; i += dimension) {
                ll = ol.proj.transform([input[i], input[i + 1]], projA, intermediate);
                point = ol.proj.transform([ll[i], ll[i + 1]], intermediate, projB);
                output[i] = point[0];
                output[i + 1] = point[1];
                for (j = dimension - 1; j >= 2; --j) {
                    output[i + j] = input[i + j];
                }
            }
        } catch (e) {}
        return output;
    };

}

if (Array.isArray(projections)) {
    for (i = 0; i < projections.length-1; i++) {
        for (j = i+1; j < projections.length; j++) {
            if (ol.proj.get(projections[i]).getCode() != ol.proj.get(projections[j]).getCode() &&
                ol.proj.get(projections[i]).getCode() != ol.proj.get(intermediate).getCode() &&
                ol.proj.get(projections[j]).getCode() != ol.proj.get(intermediate).getCode() ) {

                ol.proj.addCoordinateTransforms(
                    projections[i],
                    projections[j],
                    transform(projections[i], projections[j]),
                    transform(projections[j], projections[i])
                );

                ol.proj.addCoordinateTransforms(
                    projections[j],
                    projections[i],
                    transform(projections[j], projections[i]),
                    transform(projections[i], projections[j])
                );

            }
        }
    }
}

}

  var map = new ol.Map({
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      })
    ],
    target: 'map',
    view: new ol.View({
      projection: 'EPSG:3857',
      center: [0, 0],
      zoom: 1
    })
  });


  var queryInput = document.getElementById('epsg-query');
  var searchButton = document.getElementById('epsg-search');
  var resultSpan = document.getElementById('epsg-result');
  var renderEdgesCheckbox = document.getElementById('render-edges');

  function setProjection(code, name, proj4def, bbox) {
    if (code === null || name === null || proj4def === null || bbox === null) {
      resultSpan.innerHTML = 'Nothing usable found, using EPSG:3857...';
      map.setView(new ol.View({
        projection: 'EPSG:3857',
        center: [0, 0],
        zoom: 1
      }));
      return;
    }

    resultSpan.innerHTML = '(' + code + ') ' + name;

    var newProjCode = 'EPSG:' + code;
    if (newProjCode != 'EPSG:4269') {
      // 'EPSG:4269' is predefined in proj4 but not OL so proj4 must be registered if using OL5
      proj4.defs(newProjCode, proj4def);
    }
    if (ol.proj.proj4 && ol.proj.proj4.register) { ol.proj.proj4.register(proj4); };
    var newProj = ol.proj.get(newProjCode);
   // var fromLonLat = ol.proj.getTransform('EPSG:4326', newProj);

    if (newProjCode != 'EPSG:4326' && newProjCode != 'EPSG:3857') {
      reprojectionErrorHandler(['EPSG:4326', newProj]);
      reprojectionErrorHandler(['EPSG:3857', newProj]);
    }

    // very approximate calculation of projection extent
   // var extent = ol.extent.applyTransform(
   //     [bbox[1], bbox[2], bbox[3], bbox[0]], fromLonLat);
   // newProj.setExtent(extent);
    var newView = new ol.View({
      projection: newProj,
      center: [0, 0],
      zoom: 4
    });
    map.setView(newView);
   // newView.fit(extent);
  }


  function search(query) {
    resultSpan.innerHTML = 'Searching ...';
    fetch('https://epsg.io/?format=json&q=' + query).then(function(response) {
      return response.json();
    }).then(function(json) {
      var results = json['results'];
      if (results && results.length > 0) {
        for (var i = 0, ii = results.length; i < ii; i++) {
          var result = results[i];
          if (result) {
            var code = result['code'], name = result['name'],
                proj4def = result['proj4'], bbox = result['bbox'];
           // if (code && code.length > 0 && proj4def && proj4def.length > 0 &&
           //     bbox && bbox.length == 4) {
              setProjection(code, name, proj4def, bbox);
              return;
           // }
          }
        }
      }
      setProjection(null, null, null, null);
    });
  }


  /**
   * Handle click event.
   * @param {Event} event The event.
   */
  searchButton.onclick = function(event) {
    search(queryInput.value);
    event.preventDefault();
  };


  /**
   * Handle change event.
   */
  renderEdgesCheckbox.onchange = function() {
    map.getLayers().forEach(function(layer) {
      if (layer instanceof ol.layer.Tile) {
        var source = layer.getSource();
        if (source instanceof ol.source.TileImage) {
          source.setRenderReprojectionEdges(renderEdgesCheckbox.checked);
        }
      }
    });
  };

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

.map {
  width: 100%;
}

<link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">
<!-- The line below is only needed for old environments like Internet Explorer and Android 4.x -->
<script src="https://cdn.polyfill.io/v2/polyfill.min.js?features=requestAnimationFrame,Element.prototype.classList,URL"></script>
<script src="https://openlayers.org/en/v4.6.5/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>
<form class="form-inline">
  <label for="epsg-query">Search projection:</label>
  <input type="text" id="epsg-query" placeholder="4326, 27700, US National Atlas, Swiss, France, ..." class="form-control" size="50" />
  <button id="epsg-search" class="btn">Search</button>
  <span id="epsg-result"></span>
  <div>
    <label for="render-edges">
      Render reprojection edges
      <input type="checkbox" id="render-edges">
    </label>
  </div>
</form>

Answer 2

I was debugging why this breaks and found the first error is thrown when converting the following coord from 32661 to 3857. You can recreate this in isolation by running the following JS in a console, assuming the EPSG:32661 proj4js def is loaded:

proj4("EPSG:32661", "EPSG:3857", [2000000.0000000005, 2000000])

-> breaks

whereas, to WGS84,

proj4("EPSG:32661", "EPSG:4326", [2000000.0000000005, 2000000])

-> [90, 90]

I verified this is the same in proj.4 itself via a query in:

https://mygeodata.cloud/cs2cs/

I don't know enough to say whether this is OK or a bug. It's interesting that it happens at what is basically the north pole, so maybe there's some rounding errors or math error at this extreme point?