Filter order estimation
Theoretically it's possible to do this, although it will often not be practical.

Let's consider this in polynomial space. For a filter of order N you have 2*N+1 independent variables (N for the denominator and N+1 for the numerator). Let's look an arbitrary point $z_{k}$ in the z-plane and let's say the value of the transfer function at this point is H($z_{k}$). The relationship between the transfer function and all the filter coefficients can be written as equation that's linear in all filter coefficients as follows: 
$$\sum_{n=0}^{2*N}b_{n}\cdot z_{k}^{-n} - H(z_{k})\cdot \sum_{n=1}^{2*N}a_{n}\cdot z_{k}^{-n}=H(z_{k})$$
So if you pick M different frequencies $z_{k}$ you'll end up with a set of M complex linear equations or 2*M real equations. Since your number of unknowns is odd (2*N+1) you probably always want to pick one frequency where z is real, i.e z = 1 or $\omega$ = 0.

If M is larger than N than the system of equations is linearly dependent. You can find the filter order by starting at N=1 and increase N until the equation system becomes linearly dependent. The largest N at which the system is linearly independent is the actual filter order. For this approach, it doesn't even matter what frequencies you pick. As long as they are different, any set of frequencies will work.

However, this a numerically very tricky problem. Polynomial representation for larger filter orders is numerically very fragile and the smallest amount of noise or uncertainty leads to very large numerical errors. For example, if you determine the values of the sampled transfer function through measurement, the required measurement accuracy will be prohibitive unless it's very benign low order filter.