The numerical contour integral method with hyperbolic contour is exploited to solve space-fractional diffusion equations. By making use of the Toeplitz-like structure of spatial discretized matrices and the relevant properties, the regions that the spectra of resulting matrices lie in are derived. The resolvent norms of the resulting matrices are also shown to be bounded outside of the regions. Suitable parameters in the hyperbolic contour are selected based on these regions to solve the fractional diffusion equations. Numerical experiments are provided to demonstrate the efficiency of our contour integral methods.