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

Hi @letega. 

Thank you. Yes, certainly worth pointing out.

As I had noted previously, the code I've used here is stripped down and re-parameterized from that linked to elsewhere by user @coder1987, and the error evidently stems from there. As you have seen, record 274147 actually relates to the Three Peaks Mountain in Monterey, CA, and not Diablo. The correct entry for Mount Diablo would be record 222343.

Even then, the coordinates for this entry are given as 37.8816462, -121.9140403 and not the 37.881628, -121.914382 in the original codebase.

I have to think that, as the code was auto-generated initially using AI, the author's assistant had got muddled up a little with both its figures and its source references.

Of the various values available, I find that the coordinates given by the Wikipedia article pin the actual summit most closely, and so I have changed to these values also in the code, although the difference this makes is minimal.

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

#!/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]:
    cases: Dict[int, Tuple] = {}
    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:
    return {
        "ANCHOR_LAT": float(POINTS),
        "ANCHOR_LON": float(POINTS),
        "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:
    """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) -> Tuple:
    distance_miles = distance_inches * C
    mag_bearing_deg = (clock_hour % 12) * 30
    true_bearing_deg = (mag_bearing_deg + C) % 360
    return project_point(
        C, C,
        true_bearing_deg, distance_miles, C,
    )

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

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

if __name__ == "__main__":
    generate_conversions()

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!

MAG_DEC_1970_DEG_EAST = 12.4
MAP_SCALE_MI_PER_IN = 7.02

def generate_cases() -> Dict[int, Tuple]:
    cases: Dict[int, Tuple] = {}
    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:

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.

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.

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.

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.