2

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.

  1. Get a bunch of lat-lng coords of the US-Canada border.
  2. 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.
  3. For each of these vectors, find the unit vector perpendicular to the vector in the direction into Canada.
  4. 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:

  1. 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.

  2. 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. Any help would be greatly appreciated - many thanks in advance!

1

A simple approach would be:

  1. Load the border line into QGIS
  2. 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
  3. Remove the first layer to avoid confusion
  4. 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
  5. 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
  6. Cut the line at the westmost and eastmost point
  7. 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.

  • Andre, many thanks for this response. Since I am new to QGIS, I need more detail on how to carry out some of the steps. How would I do step #2 above? Would I right-click on the layer, select "Save as..." and then select ...? For the CRS do I select "Layer CRS" or "Project CRS"? Right one of the boxes in the save dialog box contains "WGS 84". What should I change that to? Thanks! – Shafique Jamal May 5 '14 at 12:57
  • Save As ... is correct, then Change the dropdwn box from Layer CRSto Selected CRS, which enables the second Browse button. there you can put 5070 in the filter, and then you can select NAD 83 / Conus Albers. Don't forget to choose a new file name with the first Browse button as well. And check Add saved file to map. – AndreJ May 5 '14 at 13:04
  • Hello Andre, I think I managed to get Steps 1 through 3 working now. Step 3 asked me to specify an output shapefile, (which I named latlng4). Regarding step 4, in the dialog box for Polygons to lines there is a dropdown box for Input polygon vector layer, but there are no items in this dropdown box to select, and the dialog box requires me to select something. Sorry if my questions seem overly simplistic- I'm still new to QGIS. Can you advise me on how to proceed? Thanks! – Shafique Jamal May 5 '14 at 13:31
  • In Step 3 there was a checkbox Add saved file to map which you missed. You can load the file manually (if you find it on your disk). I guess the next steps will have the same... – AndreJ May 5 '14 at 14:41
  • I think the problem I had before was that I was saving as a CSV (it was set to that from before). I changed the format to ESRI Shapefile, and now it shows up properly in the Layers sidebar. I'm using QGIS 2.2, and I changed the rendering at the bottom right corner to EPSG:5070. The coordinates at the bottom center are no longer in degrees, but look like: 1913288, 2049948. Now when I go to add a buffer, 120000 is way too high - it colors the entire canvas. A value of 1 seems more reasonable, but I don't know what distance a value of 1 corresponds to. Any thoughts on where I went wrong? – Shafique Jamal May 5 '14 at 15:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.