As the majority of us are limited in only being able to work with digital copies of the original map, we have to be careful in evaluating how to measure distances, in inches, upon it. As hinted, we are going to find that the reported scale of 6.4 inch/ml is insufficient for the task ahead, and we will be working towards obtaining a more correct value. However, there is a second related way in which it would be possible for us to go astray, if we are not careful in eliminating it at the start.

Can we be sure, in working with our scanned copies of the Phillips Map, that what the (incorrect) scale on the map claims to be an inch is, in fact, an inch?

Thankfully, we can find a comparison scan courtesy of Mike Butterfield on his website zodiackillerfacts.com that, assuming it to be precisely equivalent to the map sent by the Zodiac, allows us to check the map scale against an inch ruler at the foot of the image.

Slightly unhelpfully, although the printed scale claims 6.4 miles to the inch, the actual inch marker is not present and is, therefore, open to a degree of interpretation when trying to measure against it. To get around this, however, we can use the clearly marked 8 mile marker which, having done the sums, should correspond precisely to 1.25, or 1¼, inches.

We therefore use this marker to verify that the scale is at least accurate as to how it renders inches.

The next thing to attempt is to produce an example which is able to demonstrate the extent of the actual errors that are introduced when we fail to take into consideration the deviations, or inaccuracies, of the map we are given.

For this we will first follow the method detailed in the thread linked in the initial post and plot on the clean Phillips 66 the location corresponding to 4 inches along the 8th radian, so:

For the bearing we have heeded, as we must, the north directional marker at the bottom-left of the figure, even though this will be shown to be inaccurate. From map north we have then proceeded round through 8 × 30° + 17° = 257° to get the direction of the 8th radian, extending our range out to the equivalent of the 4 inches using the provided scale.

This takes us to somewhere around McKinley Square (perhaps just south of) off the US 101.

Next, our task is to run the appropriate math to deduce the lat./long. coordinates of that place that is:

1. on a bearing, centred on Mount Diablo, of 257° from true north (thereby ignoring, for the moment, the map’s error in orientation); and
2. at a distance from the same of 4 × 6.4 = 25.6 miles (thereby ignoring, likewise, the map’s error in scale).

For this step, however, it is probably going to be useful if we hack together a slice of code to do the work for us.

For this step, we’ll borrow a few blocks of code from the example in Python generated for D. Stampher (@coder1987) and referenced from this thread which, for our purposes, we can simplify down to:

#!/usr/bin/env python3

from typing import Dict, Tuple, Any
import math

EARTH_RADIUS_MI = 3958.8

MAG_DEC_1970_DEG_EAST = 17
MAP_SCALE_MI_PER_IN = 6.4

def generate_cases() -> Dict[int, Tuple[float, float]]:
    cases: Dict[int, Tuple[float, float]] = {}
    return {
        1: (8.0, 4.0),
    }

POINTS = {
    "mt_diablo": {
        "lat": 37.881628,
        "lon": -121.914382,
        "source": "USGS GNIS",
        "source_url": "https://edits.nationalmap.gov/apps/gaz-domestic/public/gaz-record/274127",
        "uncertainty_m": 10,
        "notes": "Anchor point for all bearings/projections."
    },
}

def build_runtime_constants() -> Dict[str, Any]:
    return {
        "ANCHOR_LAT": float(POINTS["mt_diablo"]["lat"]),
        "ANCHOR_LON": float(POINTS["mt_diablo"]["lon"]),
        "MAG_DECLINATION": float(MAG_DEC_1970_DEG_EAST),
        "MAP_SCALE": float(MAP_SCALE_MI_PER_IN),
        "EARTH_RADIUS_MI": float(EARTH_RADIUS_MI),
    }


def project_point(start_lat: float, start_lon: float,
                  bearing_deg: float, distance_miles: float,
                  earth_radius: float = 3958.8) -> Tuple[float, float]:
    """Forward geodesic projection on a sphere. Returns (lat, lon) in degrees."""
    lat1 = math.radians(start_lat)
    lon1 = math.radians(start_lon)
    brng = math.radians(bearing_deg)
    d = distance_miles / earth_radius
    lat2 = math.asin(
        math.sin(lat1) * math.cos(d) + math.cos(lat1) * math.sin(d) * math.cos(brng)
    )
    lon2 = lon1 + math.atan2(
        math.sin(brng) * math.sin(d) * math.cos(lat1),
        math.cos(d) - math.sin(lat1) * math.sin(lat2),
    )
    return math.degrees(lat2), math.degrees(lon2)

def project_from_anchor(distance_inches: float, clock_hour: int,
                        C: Dict[str, Any]) -> Tuple[float, float]:
    distance_miles = distance_inches * C["MAP_SCALE"]
    mag_bearing_deg = (clock_hour % 12) * 30
    true_bearing_deg = (mag_bearing_deg + C["MAG_DECLINATION"]) % 360
    return project_point(
        C["ANCHOR_LAT"], C["ANCHOR_LON"],
        true_bearing_deg, distance_miles, C["EARTH_RADIUS_MI"],
    )

def convert()->Dict[int, Tuple[float, float]]:
    print("Transforming map coords. to lat./long. ....\n")
    cases: Dict[int, Tuple[float, float]] = generate_cases()
    conversions: Dict[int, Tuple[float, float]] = {}
    C = build_runtime_constants()
    for case in cases:
        angle, dist = cases[case]
        conversions[case] = project_from_anchor(dist, angle, C)
    return conversions

def generate_conversions():
    conv: Dict[str, Tuple[float, float]] = convert()
    for item in conv:
        print(f"{item}: {conv[item]}")

if __name__ == "__main__":
    generate_conversions()

which covers, for now, only the single coordinate pair of 4 inches along the 8th radian.

We will be adding others later once we have made some progress with correcting the noted errors in this implementation.

When plotted onto Google Maps, the resultant lat./long. coordinates (37.79739830023115, -122.3712565315443) then bring us to the following location, which is clearly and grossly in error as, from the previous comments, we were expecting.

As a rough comparison, we can illustrate the approximate degree of the error by plotting this same location back onto the Phillips 66, so:

Again, as noted, the deviation here is seen to arise from the combined effect of neglecting both the error in orientation and that in scale, as is inherent in the Phillips Map.

We should next turn to figuring out how to correctly account for these errors in our coded model.

It is possible, of course, to treat the analysis and correction of the two distinct errors separately. We shall start with the error in orientation.

On coarse inspection, we can estimate that this amounts to the introduction of an erroneous rotation on conversion to lat./long. of around about 5° clockwise. As we have seen, this arises because the Phillips 66 is actually oriented with its Map North direction approximately this amount off of True North in the counterclockwise direction. Thus, when we apply an adjustment to Magnetic North on the map by, say, 17° E, as we are asked to do, the equivalent adjustment that we must then take to Google Maps should be not the full 17°, ignoring the error, but rather 17° – 5° = 12°.

This can be easily corrected by modifying the magnetic declination in our code accordingly, as:

MAG_DEC_1970_DEG_EAST = 12.0
MAP_SCALE_MI_PER_IN = 6.4

Running the code and plotting again to Google, we get the following as our first attempt at adjusting for the orientation error, which is actually pretty good for an initial guess by eye:

The precise error, if we can indeed determine one, is likely to be a very little less than 5°. However, we can hold that for further fine adjustments later and turn next to the error in scale.

After only a little further experimentation, with the scale in our code adjusted as:

MAG_DEC_1970_DEG_EAST = 12.0
MAP_SCALE_MI_PER_IN = 7.02

we are brought to the following, which in this one instance is now appearing pretty darn close.

Before we settle on any specific corrections, we have to check that our method of analysis is not undermined by the presence of any significant local deviations also. We will want to be confident that, having a more correct value for the scale, for example, we will find this stable right across the map, no matter where we start from or end at. Our hope is that the map is, at the very least, largely self-consistent in this sense. Still, we have to test it, and we do this as follows.

For our purposes here, we can select a small number of locations at various bearings and ranges on the map and have them each centrally mapped in the usual way from Mount Diablo. Convenience persuades us to use only whole number map distances, and we can plot our choices on the map as we have done previously.

Here, then, we will be using:

1. 5 inches along the 6th radian;
2. 2 inches along the 7th radian;
3. 3 inches along the 8th radian;
4. 4 inches along the 8th radian;
5. 5 inches along the 9th radian;
6. 3 inches along the 10th radian;
7. 1 inch along the 11th radian; and
8. 2 inches along the 12th, or 0th, radian.

Direct plotting of these eight examples on the Phillips 66 map produces the following spread.

We then move on to calculating corresponding lat./long. coordinate pairs for each of these locations, having modified our original code to include our selection, as well as some small but indicated adjustments to our corrected map parameters, as with:

MAG_DEC_1970_DEG_EAST = 12.4
MAP_SCALE_MI_PER_IN = 7.02

def generate_cases() -> Dict[int, Tuple[float, float]]:
    cases: Dict[int, Tuple[float, float]] = {}
    return {
        1: (6.0, 5.0),
        2: (7.0, 2.0),
        3: (8.0, 3.0),
        4: (8.0, 4.0),
        5: (9.0, 5.0),
        6: (10.0, 3.0),
        7: (11.0, 1.0),
        8: (12.0, 2.0),
    }

This gives us, as results:

Transforming map coords. to lat./long. ....

1: (37.38539570365206, -122.05166999719813)
2: (37.73144596298407, -122.08762901429596)
3: (37.78889318417564, -122.28202014122363)
4: (37.75772804098497, -122.40436042530442)
5: (37.98903911832634, -122.54391870173873)
6: (38.08681110021985, -122.20035444361838)
7: (37.978466303347425, -121.95335599396958)
8: (38.08007595768963, -121.85894859106457)

which points, in turn, map on Google Maps as in the following screenshot:

From a side-by-side comparison of the last two images, it’s clear that the applied corrections are very close to what is required and that there is a general consistency in how the map presents its errors.

From this we are confident in being able to state that the Phillips 66 Map, as we were given it, is inherently misaligned as to true North to an amount of approximately 4.6° W, while it is, in fact, scaled at near to 7.02 miles to the inch, not 6.4 as the map itself indicates.

And, of course, these are the corrections that must be taken into consideration, whenever we seek to convert our direct map coordinates—our latest, and bestest, killer Z32 solve—into coordinates for mapping in our modern apps.

Whether we are lining up a wooden inch ruler over the actual Phillips 66, or getting Google Maps to spit out the precise location of our latest discovery, we have to remain clear on one thing in particular—if your app solution doesn’t to some acceptable degree of accuracy match, as to location, what you would obtain from the simple ruler-on-map solution, then you’re simply not doing it right.

Because they must match.

Because, that’s how maps work!

Here’s the final code for the above test cases, in case you’d maybe like to try out a few yourself.

#!/usr/bin/env python3

from typing import Dict, Tuple, Any
import math

EARTH_RADIUS_MI = 3958.8

MAG_DEC_1970_DEG_EAST = 12.4
MAP_SCALE_MI_PER_IN = 7.02

def generate_cases() -> Dict[int, Tuple[float, float]]:
    cases: Dict[int, Tuple[float, float]] = {}
    return {
        1: (6.0, 5.0),
        2: (7.0, 2.0),
        3: (8.0, 3.0),
        4: (8.0, 4.0),
        5: (9.0, 5.0),
        6: (10.0, 3.0),
        7: (11.0, 1.0),
        8: (12.0, 2.0),
    }

POINTS = {
    "mt_diablo": {
        "lat": 37.881628,
        "lon": -121.914382,
        "source": "USGS GNIS",
        "source_url": "https://edits.nationalmap.gov/apps/gaz-domestic/public/gaz-record/274127",
        "uncertainty_m": 10,
        "notes": "Anchor point for all bearings/projections."
    },
}

def build_runtime_constants() -> Dict[str, Any]:
    return {
        "ANCHOR_LAT": float(POINTS["mt_diablo"]["lat"]),
        "ANCHOR_LON": float(POINTS["mt_diablo"]["lon"]),
        "MAG_DECLINATION": float(MAG_DEC_1970_DEG_EAST),
        "MAP_SCALE": float(MAP_SCALE_MI_PER_IN),
        "EARTH_RADIUS_MI": float(EARTH_RADIUS_MI),
    }


def project_point(start_lat: float, start_lon: float,
                  bearing_deg: float, distance_miles: float,
                  earth_radius: float = 3958.8) -> Tuple[float, float]:
    """Forward geodesic projection on a sphere. Returns (lat, lon) in degrees."""
    lat1 = math.radians(start_lat)
    lon1 = math.radians(start_lon)
    brng = math.radians(bearing_deg)
    d = distance_miles / earth_radius
    lat2 = math.asin(
        math.sin(lat1) * math.cos(d) + math.cos(lat1) * math.sin(d) * math.cos(brng)
    )
    lon2 = lon1 + math.atan2(
        math.sin(brng) * math.sin(d) * math.cos(lat1),
        math.cos(d) - math.sin(lat1) * math.sin(lat2),
    )
    return math.degrees(lat2), math.degrees(lon2)

def project_from_anchor(distance_inches: float, clock_hour: int,
                        C: Dict[str, Any]) -> Tuple[float, float]:
    distance_miles = distance_inches * C["MAP_SCALE"]
    mag_bearing_deg = (clock_hour % 12) * 30
    true_bearing_deg = (mag_bearing_deg + C["MAG_DECLINATION"]) % 360
    return project_point(
        C["ANCHOR_LAT"], C["ANCHOR_LON"],
        true_bearing_deg, distance_miles, C["EARTH_RADIUS_MI"],
    )

def convert()->Dict[int, Tuple[float, float]]:
    print("Transforming map coords. to lat./long. ....\n")
    cases: Dict[int, Tuple[float, float]] = generate_cases()
    conversions: Dict[int, Tuple[float, float]] = {}
    C = build_runtime_constants()
    for case in cases:
        angle, dist = cases[case]
        conversions[case] = project_from_anchor(dist, angle, C)
    return conversions

def generate_conversions():
    conv: Dict[str, Tuple[float, float]] = convert()
    for item in conv:
        print(f"{item}: {conv[item]}")

if __name__ == "__main__":
    generate_conversions()

And the same, here as an active and editable implementation on PlayCode, along with a few examples and modified to output direct GMaps urls:

https://playcode.io/phillips-66-to-gmaps–019e3842-51ad-76cf-9ea3-5e75bb33c715

Hi shaqmeister,

I’m not too good at maps so take this with a pinch of salt,

In the playcode link, the Source for Mt. Diablo Lat & Lon is given as: 

source_url”: “https://edits.nationalmap.gov/apps/gaz-domestic/public/gaz-record/274127”

This url gives an entry date for the data therein as “Entry Date 19 January 1981” and

Elevation: 1023 meters / 3356 feet
Coordinates: 35.8480522, -121.3054871 / 35° 50′ 52.99″ N, 121° 18′ 19.75″ W

Whereas wikipedia says[1]:
Elevation 3,849 ft (1,173 m)
Coordinates 37°52′54″N, 121°54′51″W

I don’t know how much difference that makes, but thought it worth pointing out.

1. https://en.wikipedia.org/wiki/Mount_Diablo