Exiftool will output the coordinates in low precision decimal if you use the -n switch. You can get more precision with -c switch and give the desired quantity of digits behind the decimal:
exiftool -c '%.6f' -GPSPosition filename.jpg
Shows position with 6 digits behind the decimal, which is good for finding a place within 5 inches.
Python would be so slow at this, I would use C++. Read your XYZ file using a normal ifstream then use strtok to break apart the single line into tokens and make the string tokens into double.
You will need to find the local WGS UTM coordinate reference system and project your origin to this CRS (sorry the OGR API is offline at the moment). From here you can ...
I posted a script to do this here: Exporting 3GB ArcGIS Raster to KML without losing resolution?
and since then another user posted a script here:
they both do a pretty similar thing.. reproject, tile, write a text file, package. Mine works on an image and his works on an mxd.
If you right click your layer in the table of contents, you can choose to Join data by spatial location. Alternatively, you can use the Spatial Join tool.
Do like sailors did in the past: take a mercator map that has the two locations on it, draw a straight line between the two points and measure using a protractor the angle with a meridian, that angle is your bearing.
You already found the link to maths and software
There are several projects, of various age and support status that make this faster:
Each of them accelerates spatial join, making such queries much faster, and should run them reasonably fast. You did not say ...