Thursday, 7 January 2016

latitude longitude - How do you find a lat/lon (DMS) on a paper map using a 30cm ruler?



How do you use a 30cm ruler to find DMS on a paper map? The locations I would like to find are the 'corner' points so I can generate an extent based on the four corners.


I have an old paper map (3 actually) for Northern Canada (late 1800's) that doesn't provide the Ellipsoid or the Datum. It provides a representative fraction (1:660,000 approx) and a scalebar (1" = 10 2/3 mile). The map shows grid lines spaced every 1 degree. No minutes or seconds are labelled.


I understand that NOT knowing the datum or ellipsoid will automatically introduce a margin of error into the calculations, but this is not a big deal for this exercise.


I determined the Lat/Lon of the intersecting grid lines, and from this question, was able to infer that it is closest to Lambert Conformal Conic (Statistics Canada, EPSG 3347).


Below is the index map showing all 3 maps with grid lines every 2 degrees: enter image description here


I will need to do this process for all three maps as those grid lines are spaced every 1 degree, and not 2 as in the index above.




Of course, I could geo-reference to a known spatial reference in a computerized GIS and digitize the extent, but what if your GIS is PC-less and you have traveled back in time and are now stuck...


If it is easier to provide an answer using say, an engineers ruler (1:100, 1:2500 etc) then feel free. It's just a 30cm ruler seems to be more readily available in a given situation.



Answer




This isn't so old-fashioned: I remember having to solve exactly this problem back in the 80's when we didn't have scanners readily available and had to lift coordinates and elevations off large-format printed maps for geostatistical analysis.


In effect you can already read the longitude accurately along any line of longitude on the map. You want to interpolate these measurements to four specific points (the corners). Ditto for the latitude. Thus, this problem is a special case of interpolating between contours on any contour map. Therefore you don't need to know anything about the projection or datum to do it.


Because this is supposed to be done simply, we cannot easily exploit the fact we have full contours. It will suffice to identify a few discrete points along each contour and use them. This makes the problem equivalent to the following:



Given a collection of points on map, each labeled with a (smoothly varying) numeric value, to estimate the value at one other specified point on the map.



In order to solve this we need to establish a coordinate system for the map itself. The choice doesn't matter as long as the coordinate isolines are evenly spaced (they don't even have to be mutually perpendicular!) A simple way to accomplish this is to use the ruler to measure distances from the left edge (x) and bottom edge (y) of the map. (If you have a scanned image, just use the row and column indexes of the pixels.)


The interpolation can be accomplished by fitting a trend to the data.


We know, just by looking at the map (that is, by observing the locally regular spacings of the contours), that a linear estimator will work fairly well and a quadratic estimator will work even better. It's probably overkill (and too much work) to use any higher-order estimator. A quadratic estimator requires at least six control points. Use a collection of points clustered near the estimation point: this will assure high accuracy. Use more than the minimum: this provides useful cross-checks and can even yield error estimates.


This results in the following procedure, to be done for latitude and repeated for each corner point and then repeated all over again for longitude:





  • Mark off more than six points along relevant contour lines in the vicinity of a corner point. Use several different contour levels.




  • Measure (x,y) at the marked points and at the corner point.




  • Record (x,y,dependent value) at each marked point.





  • Compute the least squares fit of the data using the model:


    (lat or lon) = a + b*x + c*y + d*x*x + e*x*y + f*y*y + error


  • Apply the fitted model to the (x,y) value for the corner point.




People have been computing least-squares fits far longer than they have had mechanical calculators available. If you really don't have a computer or calculator available, settle for a linear trend and for the (easy) calculations consult any textbook on regression published before about 1970. Otherwise, you can do the fit with a graphic calculator, a spreadsheet, or (best and easiest) any full-featured statistical package. The latter will be able to provide you a prediction interval to assess the uncertainty in the estimates.


For example, I applied this procedure twice to find (lat, lon) at the upper left corner using the marked points (red for longitude, blue for latitude, yellow for the corner):



marked map


Using obvious variable names, I obtained the predicted values with two Stata 11 commands for each calculation:


regress lat x y c.x#c.y c.x#c.x c.y#c.y if lat!=0
predict lathat
regress lon x y c.x#c.y c.x#c.x c.y#c.y if lon!=0
predict lonhat

The estimated (lat, lon) of the corner point is (61.05, -136.80). The estimated error is surprisingly large (around 0.04 degree), about twice what I would expect from the resolution of the screen image. These contour lines might not be very accurately placed.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...