How to Extract (Magnetic) XYZ Data Points from Color (Airborne Magnetic Survey) Image and Legend

legend = Import["https://i.imgur.com/onzvkyP.png"];

(* extract a slice of columns containing the colours. Rotate the legend on its side. *)
legadjusted = ImageRotate[ImageCrop[ImageTake[legend, All, {40, 45}]], -Pi/2];

(* we generate a list of points that fall in between the boxes in the legend *)
coords = Table[{x, 5}, {x, 16, 398, 9.5}];

(* illustrates where I've put the points *)
Show[legadjusted, Graphics[{Point[coords]}]]

(* extract the colours at the pixel coordinates, ignore the last few bad ones*)
colours = PixelValue[legadjusted, coords][[;; -3]];

(* create a range of values *)
values = N@Range[54443, 55845, (55845 - 54443 + 1)/Length[colours]];

(* load the image and build a NearestFunction that maps colours to indices in the colour map*)
img = Import["https://i.stack.imgur.com/bi16H.jpg"];
nf = Nearest[colours -> "Index"];

(* gets the index of the nearest colours and the corresponding value *)
col2val = Function[{c}, values[[First@nf[c]]]];


(* we are done! now get a grid of values by mapping over the image *)
datagrid = Map[col2val, ImageData[img], {2}];

(* we can even turn it into monochrome if we want *)
Image[datagrid] // ImageAdjust

(* or even make a 3D landscape *)
ListPlot3D[Reverse@datagrid, PlotTheme -> "Web"]