This function should do the trick:
GraphToFunction[img_] := Block[{words, centers, xmain, ymain, xwords, ywords, pts, f},
(*Find where all the words (ie axis numbers) are in the image and store them in a list of the form {{word,bounding rectangle}...}*)
words = TextRecognize[img, "Word", {"Text", "BoundingBox"}];
(*Replace rectangles with their center points*)
words = {ToExpression[words\[Transpose][[1]]],RegionCentroid[#] & /@ (words\[Transpose][[2]])}\[Transpose];
(*Next split words into the x and y axis labels*)
{xmain, ymain} =Commonest[#, 1] & /@ (words\[Transpose][[2]]\[Transpose]);
xwords = Select[words, #[[2, 2]] == ymain[[1]] &];ywords = Select[words, #[[2, 1]] == xmain[[1]] &];
(*Next convert the plot into a list of points using based on their \color*)
pts = PixelValuePositions[img, RGBColor["#4682B4"]];
pts = DeleteDuplicates[Mean[#] & /@ (Cases[pts, {#, _}] & /@ (pts\[Transpose][[1]]))];
(*Then rescale the points based on the values and positions of the \previosuly found words*)
pts = {(xwords[[2, 1]] - xwords[[1, 1]])/(xwords[[2, 2, 1]] -xwords[[1, 2, 1]]) (#[[1]] - xwords[[1, 2, 1]]) +xwords[[1,1]], (ywords[[2, 1]] - ywords[[1, 1]])/(ywords[[2, 2, 2]] -ywords[[1, 2, 2]]) (#[[2]] - ywords[[1, 2, 2]]) +ywords[[1, 1]]} & /@ pts;
(*Finally sort points by x value and return an interpolating \function going through all the rescaled points*)
Interpolation[SortBy[pts, #[[1]] &]]
]
Here's how it works. The first step is to automatically detect the axis labels so we can automatically scale the points later. This can be achieved using the TextRecognize[]
function in Mathematica to simultaneously find the values of the axis labels and also rectangles that bound their locations. We can then compute the centers of these rectangles so that we have a list of axis labels and their corresponding pixel positions. Next we need to figure out which labels are on the x and y axis which we can quickly do by finding the most common x and y values and grouping the points by them. Next we need to convert the graph itself into a list of points. Assuming all your graphs are of the same form (ie have the same color) we can simply use the PixelValuePositions[]
function to do this. This causes a slight issue though, because the line has width and will result in a few y values for every x value. We can get around this by simply averaging all the y values of duplicate points. Then we can take two of the points from the x-axis labels and two from the y-axis to rescale these pixel position points to their actual coordinates. And finally we can bundle up our output into a nice interpolating function so that you can simply pull any value from the graph you want.
Here is a simple example using this function to get the values at all the x-axis points (Note, in all these examples img is a variable referencing your image):
In[722]:= Evaluate[GraphToFunction[img]]/@{2000,2005,2010,2015,2020}
Out[722]= {265/44,35/8,355/88,185/88,455/528}
And here is a super cool example where we re-create the graph showing that every single point was properly found:
GraphToFunction[img];
Plot[%[x], {x, 1998, 2022}, PlotRange -> All]

And finally it's also worth noting that becuase the output is an interpolating function we can do anything we could with a standard function in Mathematica, like take its integral, derivative, etc.
In[724]:= Integrate[Evaluate[GraphToFunction[img]][x],{x,2005,2010}]//N
Out[724]= 6.1132
ResourceFunction["ExtractPlotImageData"][picture]
which extracts points from an image. $\endgroup$