1
votes

I have a mapbox example where I do the following:
1) setup a view with a pitch angle of 60 degrees
2) map the canvas corner points to get the viewing frustrum as a trapezoid
3) change the bearing by 30 degrees and get and draw the trapezoid again
4) zoom out to see the trapezoids
5) manually rotate the first trapezoid by 30 degrees
6) draw the rotated trapezoid
7) change the pitch to zero to look down at the image
When I do this, the manually rotated trapezoid (red in the image) does not match the one generated using the setBearing() call (purple in the image). It appears to be skewed improperly and I have been looking at the manual rotate code for 8 hours and cannot figure out why. Am I dealing with curvature of the earth rotated co-ordinate issues or? Can someone sort this out? Thanks!

mapboxgl.accessToken = 'pk.eyJ1IjoiZm1hY2RlZSIsImEiOiJjajJlNWMxenowNXU2MzNudmkzMndwaGI3In0.ALOYWlvpYXnlcH6sCR9MJg';

var map;

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 pointsRotate(points, cx, cy, angle){
    var radians = (Math.PI / 180) * angle;
    var cos = Math.cos(radians);
    var sin = Math.sin(radians);
    var newpoints = [];

    function rotate(x, y) {
        nx = (cos * (x - cx)) + (sin * (y - cy)) + cx,
        ny = (cos * (y - cy)) + (-sin * (x - cx)) + cy;
        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);
}

map = new mapboxgl.Map({
    container: 'map',
    style: 'mapbox://styles/mapbox/streets-v9',
    center: [-122.48610019683838, 37.82880236636284],
    zoom: 17,
    pitch: 60
});

map.on('load', function () {

  // get the center and trapezoid of the zoomed in view
  var center = map.getCenter();
  var trapezoid = getViewTrapezoid();
  
  // convert the view trapezoid to a path and add it to the view
  var trapezoid_path = convertTrapezoidToPath(trapezoid);
  addLayerToMap("viewTrapezoid",trapezoid_path,'#888',4);

  // now rotate the bearing by 30 degrees to get a second view trapezoid 
  map.setBearing(30);
  setTimeout(function() {
    var trapezoid2 = getViewTrapezoid();
    var trapezoid2_path = convertTrapezoidToPath(trapezoid2);
    addLayerToMap("viewTrapezoid2",trapezoid2_path,'#f0f',2);
    
    // return to a "top down" view and zoom out to show the trapezoids
    map.setBearing(0);
    map.setZoom(13.5);
      setTimeout(function() {
        map.flyTo({ pitch: 0 });
        
        // rotate the original view trapezoid by 30 degrees and add it to the map
        var newpath = pointsRotate(trapezoid_path,center.lng,center.lat,30);
        addLayerToMap("rotatedTrapezoid",newpath,'#f00',2);
    }, 500);
  }, 500);
  
  
 
});
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>
1
Hmmm - looks like it is a problem with longitude distances getting smaller near the poles. My calculation will only work near the equator according to this post: gis.stackexchange.com/questions/206698/…> I will try their code suggestion and update...fmacdee
Ok, so I found an algorithm for calculating the length of a longitude degree based on which degree of latitude it is at. I will give that a try... gis.stackexchange.com/questions/75528/…fmacdee

1 Answers

1
votes

Ok, so it turns out that because of longitude getting thinner as you approach the poles, you first need to convert all degree co-ordinates (lat and long) to Mercator co-ordinates that take into account the non-spherical nature of the earth. I've added two functions here that convert from lat lon to mercator and vice versa and put them into the code and the result is a trapezoid that is directly on top of the one provided using the setBearing() method. Problem solved!

mapboxgl.accessToken = 'pk.eyJ1IjoiZm1hY2RlZSIsImEiOiJjajJlNWMxenowNXU2MzNudmkzMndwaGI3In0.ALOYWlvpYXnlcH6sCR9MJg';

var map;

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 pointsRotate(points, cx, cy, angle){
    var radians = (Math.PI / 180) * angle;
    var cos = Math.cos(radians);
    var sin = Math.sin(radians);
    var newpoints = [];

    function rotate(x, y) {
        nx = (cos * (x - cx)) + (sin * (y - cy)) + cx,
        ny = (cos * (y - cy)) + (-sin * (x - cx)) + cy;
        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 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 convertDegrees2Meters(points) {
    var newpoints = [];

    for(var i=0;i<points.length;i++) {
        newpoints[i] = ll2Mercator( points[i][0], points[i][1] );
    }
    return newpoints;
}
function convertMeters2Degrees(points) {
    var newpoints = [];

    for(var i=0;i<points.length;i++) {
        newpoints[i] = Mercator2ll( points[i][0], points[i][1] );;
    }
    return newpoints;
}

map = new mapboxgl.Map({
    container: 'map',
    style: 'mapbox://styles/mapbox/streets-v9',
    center: [-122.48610019683838, 37.82880236636284],
    zoom: 17,
    pitch: 60
});

map.on('load', function () {

  // get the center and trapezoid of the zoomed in view
  var center = map.getCenter();
  var trapezoid = getViewTrapezoid();
  var center_meters = ll2Mercator(center.lng,center.lat);
  // convert the view trapezoid to a path and add it to the view
  var trapezoid_path = convertTrapezoidToPath(trapezoid);
  addLayerToMap("viewTrapezoid",trapezoid_path,'#888',4);

  // now rotate the bearing by 30 degrees to get a second view trapezoid 
  map.setBearing(30);
  setTimeout(function() {
    var trapezoid2 = getViewTrapezoid();
    var trapezoid2_path = convertTrapezoidToPath(trapezoid2);
    addLayerToMap("viewTrapezoid2",trapezoid2_path,'#f0f',2);
    
    // return to a "top down" view and zoom out to show the trapezoids
    map.setBearing(0);
    map.setZoom(13.5);
      setTimeout(function() {
        map.flyTo({ pitch: 0 });
        
        // rotate the original view trapezoid by 30 degrees and add it to the map
        var tpath_meters = convertDegrees2Meters(trapezoid_path);
        var newpath_meters = pointsRotate(tpath_meters,center_meters[0],center_meters[1],30);
        var newpath = convertMeters2Degrees(newpath_meters);
        addLayerToMap("rotatedTrapezoid",newpath,'#f00',2);
    }, 500);
  }, 500); 
});
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>