As you already have a topology for the circuit, you have several options.
One option is to hand-analyse the circuit, into an expression in jω and the values of the components. This one has particularly simple form. Start from the parallel combination of L3 and R4 in parallel, let's call that Z3. Express the impedance of L3 as jωL3. Now Z3 is L3 in parallel with R4, Z3 = 1/(1/R4 + 1/jωL3). Now add R3 in series, and put that in parallel with L2, and continue up the ladder until you finally add R1, ending up with a power series in jω. It gets tedious, but it's possible.
Now compare expression power term by power term with your input expression Z(f), where f=2πω. If you have exactly the same powers of jω in each, then match the coefficients. If they result in consistent and non-negative values for L and R, then you have yourself an exact synthesis. Unless it's an exam question that has been contrived to work exactly, it's more likely that you'll have spare power terms, or get impossible component values, in which case you're doing a best fit approximation.
An alternative is optimisation, by either repeatedly analysing your circuit in Spice for instance, and tweaking values until the response matches Z(f), or by using optimising software which does that process quickly and automatically for you.
I'd recommend doing hand optimisation with SPICE first, to explore the landscape of possible values. If you let automatic optimsation rip without understanding the regions of values in which solutions will be reasonable and so being able to constrain the optimser, then they tend to disappear off to infinity, or down local rabbit holes, without getting to solutions you're happy with.
As you have a large number of scatterplots and want to fit each of them, you obviously need an automated solution. Whether an optimisation method performs well or badly depends strongly on the details of the problem. On the one hand, your R/L ladder should be quite low order and well behaved. On the other, small values for L1 would reduce the sensitivity of the input impedance to changes in the higher numbered components, which could confound automated solving. If the scatterplots are from measurements, then the noise on the measurements could also be over-fitted and so confuse an optimser.
I'd recommend starting with something simple like the Nelder-Mead Simplex method, it's a reasonably robust non-parametric steepest descent method. You would need to construct a 'distance' function to reduce the divergence between your scatterplot and the input impedance of your trial network down to a single number. The RMS or the minimax of the absolute differences are both popular measures, though whether you base them on linear or dB units depends on where you want the best fitting to occur. Do your scatterplots have phase information? Start with an 'average' of all your scatterplots if you can define such a thing. Observe what the optimiser does as it explores the surface to see whether you're on the right track. Set upper and lower limits for all of the components if you can. When in doubt, it's better to overconstrain the optimiser and then investigate the boundaries it's pushing against, than to underconstrain it and try to figure out where it's gone and why.
Before starting with a self-directing optimisation strategy like NM, I'd recommend a more robust and simpler method such as random. If you have a plausible range for each of the values, then let the computer do the heavy lifting and throw a large selection of random values from a uniform distribution of each of the component values at an analysis. Pick the best from 1e6 runs, then do it again several times. This will give you an idea of whether the problem is well behaved, or otherwise.
