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

I've been trying to work my way through this by working through search results like:

Data extraction from a picture of a graph

Numeric data from a color map

Extract bathymetry data from a map

I can't seem to make it work.

Basically, I have a GEOTIFF of the magnetic survey flown but I don't have the raw data. The GEOTIFF is here:

And the legend is here:

What I'm trying to do is to extract the Z data from the image (ie. the value in nT). I can figure out how to get X,Y but I'm struggling to figure out how to obtain the value. Any advice would be huge. Ideally, I'd like to be able to get data from the map into a CSV I can work with. Thanks!

• Can you share the attempts you've made so far? Commented Mar 31, 2021 at 21:10
• Marco, I'm brand new to the software so I'm trying to figure it out. I get lost in the coding. I'm trying to follow this: mathematica.stackexchange.com/questions/187710/… I've imported the map and am confused with how to follow the rest: scale = Mean /@ ImageData[img][[27 ;; -26, -90 ;; -60]]; actualimage = ImageData[img][[;; , 1 ;; -112]]; nf = Nearest[scale -> Reverse@Subdivide[-7.2, 11.0, Length[scale] - 1]]; T = ArrayReshape[nf[Flatten[actualimage, 1], 1],Dimensions[actualimage][[1 ;; 2]]]; Commented Mar 31, 2021 at 21:16

1 Answer

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.sstatic.net/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"]


• +1 Why is your3d image upside down? Commented Apr 1, 2021 at 13:09
• @joojaa good point - it's because the bottom left of an image is actually the final row of the grid - (0,0) is in the top left corner when it should be at the bottom. I'll fix it. Commented Apr 1, 2021 at 15:40