One approach would be to use the frequency-domain least-squares (FDLS) method. Given a set of (complex) samples of a discrete-time system's frequency response, and a filter order chosen by the designer, the FDLS method uses linear least-squares optimization to solve for the set of coefficients (which map directly to sets of poles and zeros) for the system whose frequency response matches the desired response with minimum total squared error.
The frequency response of an $N$-th order linear discrete-time system can be written as:
$$
H(\omega) = H(z)|_{z=e^{j\omega}}
$$
where $H(z)$ is the system's transfer function in the $z$ domain. This is typically written in the rational format that follows directly from the system's difference equation:
$$
H(z) = \frac{\displaystyle\sum_{k=0}^{N}b_kz^{-k}}{1 + \displaystyle\sum_{k=1}^{N}a_kz^{-k}}
$$
The frequency response is therefore:
$$
H(\omega) = \frac{\displaystyle\sum_{k=0}^{N}b_ke^{-jk\omega}}{1 + \displaystyle\sum_{k=1}^{N}a_ke^{-jk\omega}}
$$
Rearrange the above to get:
$$
\sum_{k=0}^{N}b_ke^{-jk\omega} - H(\omega) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega}\right) = 0
$$
This equation is linear in the $2N+1$ unknown system difference equation coefficients $b_k$ and $a_k$. Given a desired frequency response $H(\omega)$, we would like to find coefficients that meet the above equation exactly for all values of $\omega$. For the general case, that's hard. So instead, we will search for a set of coefficients for a system whose frequency response approximates the desired response at a discrete set of frequencies.
In order to solve for an appropriate set of coefficients using the linear least squares method, we generate an overdetermined system of equations in those unknowns. To generate those equations, choose a collection of frequencies $\omega_m \in [0, 2\pi), m = 0, 1, \ldots , M-1$ (where $M > 2N+1$, and often $M \gg 2N+1)$. For each frequency, substitute the corresponding value of $\omega_k$ into the above equation to yield:
$$
\sum_{k=0}^{N}b_ke^{-jk\omega_k} - H(\omega_k) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega_k}\right) = 0
$$
The values $H(\omega_k)$ are obtained by sampling the desired frequency response at the chosen frequencies $\omega_k$. After generating the system of linear equations, the least-squares solution for the system's coefficients $b_k$ and $a_k$ (and therefore its poles and zeros) is easily obtained. If you substitute those coefficients back into the equation for $H(\omega)$ shown above, it should (hopefully) yield a function that is close to the template frequency response that you started with.
This technique has a few advantages:
Any arbitrary complex (magnitude and phase) frequency response can be used as the template. If you only have a magnitude constraint, you could just pick a phase response, such as linear phase.
It can be used to design both FIR and IIR filters; for an FIR realization, just remove the $a_k$ coefficients from the above.
The technique is very simple to implement and is easily parametrizable based upon the desired system order.
While there may not be a good way to estimate a priori what the required system order is to meet your design constraints, it is simple to iteratively increase the order $N$ until some selected error metric is met (such as peak error, total squared error, or deviation within a specific band).
You could extend this method a bit to use weighted least-squares optimization if needed; this would allow you to specify regions of the frequency response whose approximation error is weighted more than others. This allows you to more tightly control passband/stopband areas while allowing more slop in "don't-care" areas.
