[GIS] Plotting line on map, or obtaining latitude and longitude coordinates of such line, that is 120km north of US-Canada border

How do I plot a line on a map, or obtain the latitude and longitude coordinates of such a line, that is 120km north of the US-Canada border, including parts of the border that cross water?

I managed to get coordinates along the land border, but not along the parts of the border that cross water. I have been trying to figure this one out for a couple of days. Here are just some of the many posts I have consulted:

On how to get the lat & lng coordinates along the border:

https://stackoverflow.com/questions/12194046/google-maps-api-v3-how-to-get-region-border-coordinates-polyline-data

https://stackoverflow.com/questions/9706484/add-search-area-outline-onto-google-maps-result

These ones look relevant from the title, but I was unable to figure out how to use them… and I'm not exactly looking to create a parallel line; the new line would not be exactly parallel to the original one, though it would look so in many parts:

How to create a line in a given distance to an existing one?

Given a flight path given in latitudes and longitudes, how do I find another flight path that is parallel to the path given?

Here is the approach that I have tried.

Get a bunch of lat-lng coords of the US-Canada border.
For each pair of adjacent points, calculate a vector describing a line that joins the two points, and calculate the lat-lng coords of the midpoint of this line.
For each of these vectors, find the unit vector perpendicular to the vector in the direction into Canada.
For each unit vector, calculate the lat-lng coords of the point 120 km into Canada from the midpoint calculated earlier by adding the unit vector times distance of 120 km to the midpoint.

Regarding step 1: I got border data from the National Geographic Data Center. I loaded the appropriate .shp file into qgis. Using the QChainage function (which I found from here) I got the lat-lng coordinates of the border (mostly along the land, not on the shores near Vancouver), and then saved them as a CSV file.

For 2, 3, and 4, I tried my best, as I hope you can see from my code below:

function initialize() {  
    var locations1 = [
        [61.48674356, -141.0012206],
        [61.58674194, -141.001389],
        [61.68674007, -141.0012813],
        [61.78673774, -141.0012812],
        [61.88672508, -141.0010658],
        [61.98671795, -141.0012186],
        [62.08670822, -141.001111],
        [62.18669614, -141.0007711],
        [62.28668686, -141.0007712],
        [62.38667513, -141.0003868],
        [62.4866706, -141.0008994],
        [62.58666236, -141.0004485],
        [62.68666005, -141.0004485],
        [62.78664608, -141.0012383],
        [62.88663255, -141.0017758],
        [62.9866291, -141.0014996],
        [63.08662212, -141.0014996],
        [63.18661862, -141.0015582],
        [63.2866163, -141.0015568],
        [63.38661351, -141.001667],
        [63.48661052, -141.0014997],
        [63.58660488, -141.0017284],
        [63.68660008, -141.0015566],
        [63.78659307, -141.0017752],
        [63.88658963, -141.0014999],
        [63.98657671, -141.0004492],
        [64.08656179, -141.0007687],
        [64.18655702, -141.0003869],
        [64.28655251, -141.0008974],
        [64.38653864, -141.0007691],
        [64.48652245, -141.000774],
        [64.58650502, -141.0010629],
        [64.6864979, -141.0012204],
        [64.78649628, -141.001389],
        [64.88649628, -141.001389],
        [64.98649325, -141.0015581],
        [65.08648629, -141.0017761],
        [65.18644235, -141.000449],
        [65.2864275, -141.0007737],
        [65.38641012, -141.0010605],
        [65.48640297, -141.0012212],
        [65.58639948, -141.0014983],
        [65.6863939, -141.0017249],
        [65.78634859, -141.0003877],
        [65.88633578, -141.0012194],
        [65.98633416, -141.001389],
        [66.08632904, -141.0017238],
        [66.18630654, -141.0003881],
        [66.28629372, -141.0012198],
        [66.38629025, -141.0014978],
        [66.48628678, -141.0015563],
        [66.58626503, -141.0007817],
        [66.6862579, -141.0012781],
        [66.78625233, -141.0017225],
        [48.09135328, -89.08774781],
        [48.12990592, -88.99551236],
        [48.16846578, -88.90327996],
        [48.20654583, -88.81083824],
        [48.24454325, -88.71836407],
        [48.26866166, -88.62203098],
        [48.28553203, -88.52349756],
        [48.30237648, -88.42497186],
        [48.29773516, -88.32843646],
        [48.26155109, -88.23525283],
        [48.22619115, -88.14176979],
        [48.19037419, -88.04846379],
        [48.15399495, -87.95536856],
        [48.11835363, -87.86199437],
        [48.08199413, -87.76889477],
        [48.04571355, -87.67578091],
        [48.00961195, -87.5825905],
        [47.9730714, -87.48952847],
        [47.93657073, -87.3964743],
        [47.90027541, -87.30336795],
        [47.86346579, -87.21042154],
        [47.8267048, -87.11745782],
        [47.79009911, -87.02444783],
        [47.75299358, -86.93163457],
        [47.7161161, -86.83874347],
        [47.67919225, -86.74585927],
        [47.64177662, -86.65316338],
        [47.60480505, -86.56032087],
        [47.56753398, -86.46757675],
        [47.52968697, -86.37503409],
        [47.4923828, -86.28233453],
        [47.45496578, -86.18964134],
        [47.416985, -86.09720241],
        [47.37920703, -86.00469466],
        [47.34180458, -85.91203491],
        [47.30337016, -85.81978872],
        [47.26558586, -85.72725087],
        [47.22746189, -85.63489552],
        [47.1891533, -85.54258762],
        [47.15107131, -85.45017422],
        [47.11269326, -85.35793446],
        [47.0740636, -85.26578074],
        [47.03580892, -85.17347346],
        [46.99727296, -85.08125112],
        [46.95816275, -84.98928576],
        [46.91957328, -84.89712782],
        [46.84871995, -84.84002049],
        [46.75400677, -84.80847183],
        [46.65922855, -84.77728631],
        [46.59130709, -84.70495143],
        [46.5248768, -84.63023577],
        [46.4656176, -84.55266116],
    ];
    // I deleted a lot of points above to save space here in the post

// This is for testing
locations2 = [
            [43.7215978, -79.705822]  
];

centerpoint1 = new google.maps.LatLng(46.4656176, -84.55266116);
centerpoint2 = new google.maps.LatLng(43.7215978, -79.705822);

var x = 1;
if (x==1) {
  var centerpoint = centerpoint1;
  var locations = locations1;
} else {
  var centerpoint = centerpoint2;
  var locations = locations2;
};

var mapOptions = {
    zoom: 10,
    center: centerpoint,
    mapTypeId: google.maps.MapTypeId.ROADMAP,
    scaleControl: true
};

var map = new google.maps.Map(document.getElementById('map-canvas'), mapOptions);
var marker, i;

for (i = 0; i < locations.length; i++) {  
    marker = new google.maps.Marker({
    position: new google.maps.LatLng(locations[i][0], locations[i][3]),
    map: map
  });
}

// **********************************************************************************************************************
// Now plot markers X km from the above points
// **********************************************************************************************************************


var R = 6371; // km

function toRadians(Value) {
    /** Converts numeric degrees to radians */
    return Value * Math.PI / 180;
}

function fromRadians(Value) {
    return Value * 180 / Math.PI;
}

// find a perpendicular vector between pairs of points to get to X km from border
for (i = 0; i < locations.length - 1; i++) {
    // v1 should be the unit vector between adjacent points, and v2 will be the unit vector perpendicular to this in the direction of farther into Canada
    var v1 = [locations[i+1][0] - locations[i][0],locations[i+1][4] - locations[i][5]];

// From http://www.movable-type.co.uk/scripts/latlong.html
var lat1 = locations[i][0];
var lat2 = locations[i+1][0];
var lon1 = locations[i][6];
var lon2 = locations[i+1][7];

var lambda_1 = toRadians(lon1);

var phi_1 = toRadians(lat1);
var phi_2 = toRadians(lat2);
var delta_phi = toRadians((lat2-lat1));
var delta_lambda = toRadians((lon2-lon1));

var Bx = Math.cos(phi_2) * Math.cos(delta_lambda);
var By = Math.cos(phi_2) * Math.sin(delta_lambda);
var phi_3 = Math.atan2(Math.sin(phi_1) + Math.sin(phi_2),
                Math.sqrt( (Math.cos(phi_1)+Bx)*(Math.cos(phi_1)+Bx) + By*By ) );
var lambda_3 = lambda_1 + Math.atan2(By, Math.cos(phi_1) + Bx);

// Get the midpoint between 
midpoint_p1_p2 = [fromRadians(phi_3), fromRadians(lambda_3)];

// find a unit vector perpendicular to the line formed by the two points. Am I doing this right? Most certainly not!
var d = 1.2; // d should correspond to 120 km... still need to figure out how to do this
var y2 = lon2-lon1;
var x2 = lat2-lat1;
var k = (y2*y2)/(x2*x2);
var x1 = Math.sqrt(d*k/(1+k));
var y1 = Math.sqrt(Math.abs(d*d-x1*x1));
var v1 = [x1, y1];
var P3 = [midpoint_p1_p2[0]+x1, midpoint_p1_p2[1]+y1];

marker = new google.maps.Marker({
    position: new google.maps.LatLng(midpoint_p1_p2[0], midpoint_p1_p2[1]),
    icon: {
      path: google.maps.SymbolPath.CIRCLE,
      scale: 2
    },
    map: map
});


marker = new google.maps.Marker({
    position: new google.maps.LatLng(P3[0], P3[1]),
    icon: {
      path: google.maps.SymbolPath.CIRCLE,
      scale: 4
    },
    map: map
});

console.log('midpoint');
console.log(midpoint_p1_p2);
console.log('P3');
console.log(P3);
console.log('v1');
console.log(v1);
console.log('k');
console.log(k);

}

}

google.maps.event.addDomListener(window, 'load', initialize);

Well, this produces something:

Map with border points, midpoints and attempted calculation of points 120 km into Canada

Note: the small circles are the midpoints, the big circles are supposed to be the points going 120 km into Canada, and the red icons are the points along the border.

Problems:

I have to make sure that the calculated points 120 km into Canada are going in the right direction (into Canada) and the right distance (120 km). Sometimes the points go in the wrong direction, and I don't know how to get the distance right.
Because of the vertices in the border, the calculated 120km points will form lines that cross. I would need to figure out an algorithm for going through the points to eliminate some to get rid of this.

I can program in Python and Javascript, I know how to use the Google Maps v3 API, and I'm just starting out with QGIS.

Best Answer

A simple approach would be:

Load the border line into QGIS
Convert to a projection with meters as units using Rightclick -> Save As ..., choose shapefile format, file name, a suitable CRS (see below) and check Add saved file to map
Remove the first layer to avoid confusion
Create a buffer with a value of 120000 with Vector -> Geoprocessing Tools -> Buffer(s), choose another file name and check Add saved file to map again
Convert the resulting polygon to lines with Vector -> Geometry Tools -> Polygons to lines, choose another file name and check Add saved file to map again
Cut the line at the westmost and eastmost point
Delete the southern part

The projection should span the whole North America, so I suggest EPSG:5070 as projection. Google Mercator (or any other Mercator projection) is not suitable, as the units are not real metres. UTM has only a limited range of longitudes.

Best Answer

Related Solutions

[GIS] Accuracy of Google Maps latitude and longitude

[GIS] Creating set of equally spaced points along line using QChainage and exporting these as latitude and longitude coordinates in degrees

Related Question