Practical Least Squares FIR filter design: How to find a realistic desired response
In general you seem to be doing all the right things.

The linear phase term is required since an arbitrary magnitude and phase response is typically non-causal. That means it's impossible to fit with a causal filter and your fitting error is very high. 

One way to estimate a "good" delay value is to sample your target on an FFT grid that's relatively dense and do an inverse FFT. Ideally this looks like an actual impulse response with a clear peak someplace and trailing off at the edge. Circulate this until you've captured the bulk of the energy. 

The other options is to try brute force: just try a lot of delays and pick the one that gives you the smallest error. 

Controlling the out of band behavior you can play around with the weights of the 
"don't care" bands. It's probably best to reduce the weights until the out of band energy is just about the same as the in-band energy. If you want real bandpass behavior, just apply a regular bandpass on top of your least square design. Trying to get a decent high pass at 250 Hz will cost you a lot of coefficients otherwise.

Your target has a lot of detail at low frequencies which determines the minimum number of FIR coefficients. You could try to prefit the major low frequency wiggles with an IIR filter and then just use an FIR for the remaining difference to the target.