mapbox-gl-js: Adjust visible area & bearing to a given line, for a given pitch

Question

I'm trying to optimize a Mapbox view for long-distance hiking trails, like the Appalachian Trail or the Pacific Crest Trail. Here's an example, which I've oriented by hand, showing the Senda Pirenáica in Spain:

The area of interest, the viewport, and the pitch are given. I need to find the correct center, bearing, and zoom.

The map.fitBounds method doesn't help me here because it assumes pitch=0 and bearing=0.

I've done some poking around and this seems to be a variation of the smallest surrounding rectangle problem, but I'm stuck on a couple of additional complications:

1. How do I account for the distorting effect of pitch?
2. How do I optimize for the aspect ratio of the viewport? Note that taking the viewport narrower or wider would change the bearing of the best solution:

FWIW I'm also using turf-js, which helps me get the convex hull for the line.

Answer 1

This solution results in the path displayed at the correct bearing with a magenta trapezoid outline showing the target "tightest trapezoid" to show the results of the calculations. The extra line coming from the top corner shows where the map.center() value is located.

The approach is as follows:

1. render the path to the map using the "fitbounds" technique to get an approximate zoom level for the "north up and pitch=0" situation
   
2. rotate the pitch to the desired angle
   
3. grab the trapezoid from the canvas
   

This result would look like this:

After this, we want to rotate that trapezoid around the path and find the tightest fit of the trapezoid to the points. In order to test for the tightest fit it is easier to rotate the path rather than the trapezoid so I have taken that approach here. I haven't implemented a "convex hull" on the path to minimize the number of points to rotate but that is something that can be added as an optimization step.
To get the tightest fit, the first step is to move the map.center() so that the path is at the "back" of the view. This is where the most space is in the frustum so it will be easy to manipulate it there:

Next, we measure the distance between the angled trapezoid walls and each point in the path, saving the closest points on both the left and right sides. We then center the path in the view by translating the view horizontally based on these distances, and then scale the view to eliminate that space on both sides as shown by the green trapezoid below:

The scale used to get this "tightest fit" gives us our ranking for whether this is the best view of the path. However, this view may not be the best visually since we pushed the path to the back of the view to determine the ranking. Instead, we now adjust the view to place the path in the vertical center of the view, and scale the view triangle larger accordingly. This gives us the magenta colored "final" view desired:

Finally, this process is done for every degree and the minimum scale value determines the winning bearing, and we take the associated scale and center position from there.

mapboxgl.accessToken = 'pk.eyJ1IjoiZm1hY2RlZSIsImEiOiJjajJlNWMxenowNXU2MzNudmkzMndwaGI3In0.ALOYWlvpYXnlcH6sCR9MJg';    var map;    var myPath = [          [-122.48369693756104, 37.83381888486939],          [-122.48348236083984, 37.83317489144141],          [-122.48339653015138, 37.83270036637107],          [-122.48356819152832, 37.832056363179625],          [-122.48404026031496, 37.83114119107971],          [-122.48404026031496, 37.83049717427869],          [-122.48348236083984, 37.829920943955045],          [-122.48356819152832, 37.82954808664175],          [-122.48507022857666, 37.82944639795659],          [-122.48610019683838, 37.82880236636284],          [-122.48695850372314, 37.82931081282506],          [-122.48700141906738, 37.83080223556934],          [-122.48751640319824, 37.83168351665737],          [-122.48803138732912, 37.832158048267786],          [-122.48888969421387, 37.83297152392784],          [-122.48987674713133, 37.83263257682617],          [-122.49043464660643, 37.832937629287755],          [-122.49125003814696, 37.832429207817725],          [-122.49163627624512, 37.832564787218985],          [-122.49223709106445, 37.83337825839438],          [-122.49378204345702, 37.83368330777276]      ];    var myPath2 = [          [-122.48369693756104, 37.83381888486939],          [-122.49378204345702, 37.83368330777276]      ];    function addLayerToMap(name, points, color, width) {      map.addLayer({          "id": name,          "type": "line",          "source": {              "type": "geojson",              "data": {                  "type": "Feature",                  "properties": {},                  "geometry": {                      "type": "LineString",                      "coordinates": points                  }              }          },          "layout": {              "line-join": "round",              "line-cap": "round"          },          "paint": {              "line-color": color,              "line-width": width          }      });  }  function Mercator2ll(mercX, mercY) {       var rMajor = 6378137; //Equatorial Radius, WGS84      var shift  = Math.PI * rMajor;      var lon    = mercX / shift * 180.0;      var lat    = mercY / shift * 180.0;      lat = 180 / Math.PI * (2 * Math.atan(Math.exp(lat * Math.PI / 180.0)) - Math.PI / 2.0);        return [ lon, lat ];  }    function ll2Mercator(lon, lat) {      var rMajor = 6378137; //Equatorial Radius, WGS84      var shift  = Math.PI * rMajor;      var x      = lon * shift / 180;      var y      = Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);      y = y * shift / 180;        return [ x, y ];  }    function convertLL2Mercator(points) {      var m_points = [];      for(var i=0;i<points.length;i++) {          m_points[i] = ll2Mercator( points[i][0], points[i][1] );      }      return m_points;  }  function convertMercator2LL(m_points) {      var points = [];      for(var i=0;i<m_points.length;i++) {          points[i] = Mercator2ll( m_points[i][0], m_points[i][1] );;      }      return points;  }  function pointsTranslate(points,xoff,yoff) {      var newpoints = [];      for(var i=0;i<points.length;i++) {          newpoints[i] = [ points[i][0] + xoff, points[i][1] + yoff ];      }      return(newpoints);  }    // note [0] elements are lng [1] are lat  function getBoundingBox(arr) {      var ne = [ arr[0][0] , arr[0][1] ];       var sw = [ arr[0][0] , arr[0][1] ];       for(var i=1;i<arr.length;i++) {          if(ne[0] < arr[i][0]) ne[0] = arr[i][0];          if(ne[1] < arr[i][1]) ne[1] = arr[i][1];          if(sw[0] > arr[i][0]) sw[0] = arr[i][0];          if(sw[1] > arr[i][1]) sw[1] = arr[i][1];      }      return( [ sw, ne ] );  }    function pointsRotate(points, cx, cy, angle){      var radians = angle * Math.PI / 180.0;      var cos = Math.cos(radians);      var sin = Math.sin(radians);      var newpoints = [];        function rotate(x, y) {          var nx = cx + (cos * (x - cx)) + (-sin * (y - cy));          var ny = cy + (cos * (y - cy)) + (sin * (x - cx));          return [nx, ny];      }      for(var i=0;i<points.length;i++) {          newpoints[i] = rotate(points[i][0],points[i][1]);      }      return(newpoints);  }    function convertTrapezoidToPath(trap) {      return([           [trap.Tl.lng, trap.Tl.lat], [trap.Tr.lng, trap.Tr.lat],           [trap.Br.lng, trap.Br.lat], [trap.Bl.lng, trap.Bl.lat],           [trap.Tl.lng, trap.Tl.lat] ]);  }    function getViewTrapezoid() {      var canvas = map.getCanvas();      var trap = {};        trap.Tl = map.unproject([0,0]);      trap.Tr = map.unproject([canvas.offsetWidth,0]);      trap.Br = map.unproject([canvas.offsetWidth,canvas.offsetHeight]);      trap.Bl = map.unproject([0,canvas.offsetHeight]);        return(trap);  }    function pointsScale(points,cx,cy, scale) {      var newpoints = []        for(var i=0;i<points.length;i++) {          newpoints[i] = [ cx + (points[i][0]-cx)*scale, cy + (points[i][1]-cy)*scale ];      }      return(newpoints);  }    var id = 1000;  function convertMercator2LLAndDraw(m_points, color, thickness) {      var newpoints = convertMercator2LL(m_points);      addLayerToMap("id"+id++, newpoints, color, thickness);  }    function pointsInTrapezoid(points,yt,yb,xtl,xtr,xbl,xbr) {      var str = "";      var xleft = xtr;      var xright = xtl;        var yh = yt-yb;      var sloperight = (xtr-xbr)/yh;      var slopeleft = (xbl-xtl)/yh;        var flag = true;        var leftdiff = xtr - xtl;      var rightdiff = xtl - xtr;        var tmp = [ [xtl, yt], [xtr, yt], [xbr,yb], [xbl,yb], [xtl,yt] ];  //    convertMercator2LLAndDraw(tmp, '#ff0', 2);        function pointInTrapezoid(x,y) {          var xsloperight = xbr + sloperight * (y-yb);          var xslopeleft = xbl - slopeleft * (y-yb);            if((x - xsloperight) > rightdiff) {              rightdiff = x - xsloperight;              xright = x;          }          if((x - xslopeleft) < leftdiff) {              leftdiff = x - xslopeleft;              xleft = x;          }            if( (y<yb) || (y > yt) ) {              console.log("y issue");          }          else if(xsloperight < x) {              console.log("sloperight");          }          else if(xslopeleft > x) {              console.log("slopeleft");          }           else return(true);          return(false);      }        for(var i=0;i<points.length;i++) {          if(pointInTrapezoid(points[i][0],points[i][1])) {              str += "1";          }          else {              str += "0";              flag = false;          }      }      if(flag == false) console.log(str);        return({ leftdiff: leftdiff, rightdiff: rightdiff });  }    var viewcnt = 0;  function calculateView(trap, points, center) {      var bbox = getBoundingBox(points);      var bbox_height = Math.abs(bbox[0][1] - bbox[1][1]);      var view = {};        // move the view trapezoid so the path is at the far edge of the view      var viewTop = trap[0][1];      var pointsTop = bbox[1][1];      var yoff = -(viewTop - pointsTop);         var extents = pointsInTrapezoid(points,trap[0][1]+yoff,trap[3][1]+yoff,trap[0][0],trap[1][0],trap[3][0],trap[2][0]);        // center the view trapezoid horizontally around the path      var mid = (extents.leftdiff - extents.rightdiff) / 2;        var trap2 = pointsTranslate(trap,extents.leftdiff-mid,yoff);        view.cx = trap2[5][0];      view.cy = trap2[5][1];        var w = trap[1][0] - trap[0][0];      var h = trap[1][1] - trap[3][1];        // calculate the scale to fit the trapezoid to the path      view.scale = (w-mid*2)/w;        if(bbox_height > h*view.scale) {          // if the path is taller than the trapezoid then we need to make it larger          view.scale = bbox_height / h;      }      view.ranking = view.scale;        var trap3 = pointsScale(trap2,(trap2[0][0]+trap2[1][0])/2,trap2[0][1],view.scale);        w = trap3[1][0] - trap3[0][0];      h = trap3[1][1] - trap3[3][1];      view.cx = trap3[5][0];      view.cy = trap3[5][1];        // if the path is not as tall as the view then we should center it vertically for the best looking result      // this involves both a scale and a translate      if(h > bbox_height) {          var space = h - bbox_height;          var scale_mul = (h+space)/h;          view.scale = scale_mul * view.scale;          cy_offset = space/2;                        trap3 = pointsScale(trap3,view.cx,view.cy,scale_mul);                trap3 = pointsTranslate(trap3,0,cy_offset);          view.cy = trap3[5][1];      }        return(view);  }    function thenCalculateOptimalView(path) {      var center = map.getCenter();      var trapezoid = getViewTrapezoid();      var trapezoid_path = convertTrapezoidToPath(trapezoid);      trapezoid_path[5] = [center.lng, center.lat];        var view = {};      //addLayerToMap("start", trapezoid_path, '#00F', 2);        // get the mercator versions of the points so that we can use them for rotations      var m_center = ll2Mercator(center.lng,center.lat);      var m_path = convertLL2Mercator(path);      var m_trapezoid_path = convertLL2Mercator(trapezoid_path);        // try all angles to see which fits best      for(var angle=0;angle<360;angle+=1) {          var m_newpoints = pointsRotate(m_path, m_center[0], m_center[1], angle);          var thisview = calculateView(m_trapezoid_path, m_newpoints, m_center);          if(!view.hasOwnProperty('ranking') || (view.ranking > thisview.ranking)) {                         view.scale = thisview.scale;              view.cx = thisview.cx;              view.cy = thisview.cy;              view.angle = angle;              view.ranking = thisview.ranking;          }      }        // need the distance for the (cx, cy) from the current north up position      var cx_offset = view.cx - m_center[0];       var cy_offset = view.cy - m_center[1];      var rotated_offset =  pointsRotate([[cx_offset,cy_offset]],0,0,-view.angle);        map.flyTo({ bearing: view.angle, speed:0.00001 });        // once bearing is set, adjust to tightest fit      waitForMapMoveCompletion(function () {          var center2 = map.getCenter();          var m_center2 = ll2Mercator(center2.lng,center2.lat);          m_center2[0] += rotated_offset[0][0];                  m_center2[1] += rotated_offset[0][1];          var ll_center2 = Mercator2ll(m_center2[0],m_center2[1]);          map.easeTo({              center:[ll_center2[0],ll_center2[1]],               zoom : map.getZoom() });          console.log("bearing:"+view.angle+ " scale:"+view.scale+" center: ("+ll_center2[0]+","+ll_center2[1]+")");            // draw the tight fitting trapezoid for reference purposes              var m_trapR = pointsRotate(m_trapezoid_path,m_center[0],m_center[1],-view.angle);          var m_trapRS = pointsScale(m_trapR,m_center[0],m_center[1],view.scale);          var m_trapRST = pointsTranslate(m_trapRS,m_center2[0]-m_center[0],m_center2[1]-m_center[1]);          convertMercator2LLAndDraw(m_trapRST,'#f0f',4);      });  }    function waitForMapMoveCompletion(func) {      if(map.isMoving())           setTimeout(function() { waitForMapMoveCompletion(func); },250);      else          func();  }    function thenSetPitch(path,pitch) {      map.flyTo({ pitch:pitch } );      waitForMapMoveCompletion(function() { thenCalculateOptimalView(path); })  }    function displayFittedView(path,pitch) {      var bbox = getBoundingBox(path);      var path_cx = (bbox[0][0]+bbox[1][0])/2;      var path_cy = (bbox[0][1]+bbox[1][1])/2;        // start with a 'north up' view      map = new mapboxgl.Map({          container: 'map',          style: 'mapbox://styles/mapbox/streets-v9',          center: [path_cx, path_cy],          zoom: 12      });        // use the bounding box to get into the right zoom range      map.on('load', function () {          addLayerToMap("path",path,'#888',8);          map.fitBounds(bbox);          waitForMapMoveCompletion(function() { thenSetPitch(path,pitch); });      });  }    window.onload = function(e) {      displayFittedView(myPath,60);  }

body { margin:0; padding:0; }  #map { position:absolute; top:0; bottom:0; width:100%; }

<script src='https://api.tiles.mapbox.com/mapbox-gl-js/v0.37.0/mapbox-gl.js'></script>  <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.37.0/mapbox-gl.css' rel='stylesheet' />  <div id='map'></div>