Bilinear transform with pre-warping for systems other than classical filters
The bilinear transform ("Tustin's method") is indeed mainly used for transforming frequency selective filters with magnitude responses that are optimal with respect to some criterium, such as Butterworth, Chebyshev, or Cauer filters.

For more general systems, the bilinear transform is usually not the best choice because of the frequency warping, which mainly affects frequency domain behavior close to Nyquist. For frequency-selective filters, cut-off frequencies can be pre-warped such that the cut-off frequencies of the resulting discrete-time system are correct, but for more general systems it is usually not clear which frequencies should be matched and which shouldn't. Also, for filters with approximately piece-wise constant magnitude responses, frequency warping is not a big problem as long as the cut-off frequencies are mapped correctly. For general systems one would usually prefer the error between the continuous-time system and the discrete-time system to be distributed more evenly in the frequency domain, which is not possible with the bilinear transform.

There are several methods for converting continuous-time systems to discrete-time systems. All of them have advantages and disadvantages, and it depends on the application which one is to be preferred. Take a look at Matlab's c2d function, which implements several methods. Also take a look at this question and its answers.

For approximating a measured response, you would directly design the corresponding system in the discrete-time domain. For general frequency responses, the frequency sampling method is very simple and effective. If an IIR approximation is desired, the Matlab function invfreqz could be used. It tries to fit a rational transfer function to given frequency response data.