3

I am trying to fit some horizontal wind data to a cosine curve in order to estimate the direction and speed of winds at various heights (Velocity Azimuth Display), however it appears that whenever I attempt to do so with values > ~1, the curve appears way too flat and the output of the fit is lower than expected. 

import numpy as np
import scipy.optimize as sc

azimuth = np.full((8), 60) #All values = 60 deg.
velocity = [5.6261001,6.6962662,3.9316666,-0.88413334,-5.4323335,-6.5153003,-3.2538002,1.0269333]
#Function that defines curve that data will be fitted to
def cos_Wave(x,a, b, c):
    return a * np.cos(x-b) + c

azimuthData = np.deg2rad(azimuth)
coeffs, matcov = sc.curve_fit(cos_Wave, azimuthData, velocity, p0 = (1,0,0)

plt.scatter(azimuthData, velocity)
plt.plot(azimuthData, cos_Wave(azimuthData, *coeffs))
plt.show()
print(coeffs)

With the coeffs output being:[ 1., 0., 0.13705066] and plot attached:

Python CurveFit

I have performed a similar curvefit using IDL's builtin curvefit function, and received more realistic values yielding [ 7.0348234, 0.59962606, 0.079354301] and providing a proper fit. Is there any reason why this is the case? I am assuming that it likely has something to do with the initial estimate (P0), however, utilizing an initial initial estimate in the IDL implementation still provides much more reasonable results.

asynchronos
  • 583
  • 2
  • 14
bgoudeau
  • 33
  • 4
  • `azimuth = [0:8]` is not valid python. Maybe you're trying to simplify the code for the question, but that is not a good way to do it. Ideally, your code is something we can copy and run ourselves. – Warren Weckesser Jun 26 '17 at 22:44
  • Sorry about that, corrected – bgoudeau Jun 27 '17 at 00:12
  • There is more to be fixed. The line containing the call to `curve_fit` is missing a closing parenthesis. More importantly, `azimuth` (and therefore `azimuthData`) is an array containing all the same values. If these are the *x* coordinates of your data, then of course `curve_fit` won't work. However, I think this code is not the code that generated the plot that you show. Please fix the code, and verify that the code that you include in the question is exactly the code that you have a question about. That is, run *this code*, and make sure it demonstrate the problem that you have. – Warren Weckesser Jun 27 '17 at 02:23
  • @bgoudeau All your azimuth values are the same. Of course, you will get a linear fit instead of a wave. Then, the `p0` parameter in `curve_fit` takes in a list rather than a tuple. – David Jaimes Jun 27 '17 at 05:53

1 Answers1

3

You need to fix a few things:

import numpy as np
import scipy.optimize as sc
import matplotlib.pyplot as plt

azimuth = np.linspace(0, 360, 8)  # deg values from 0 to 360
velocity = [5.6261001, 6.6962662, 3.9316666, -0.88413334, -5.4323335,
            -6.5153003, -3.2538002, 1.0269333]


def cos_Wave(x, a, b, c):
    """Function that defines curve that data will be fitted to"""
    return a * np.cos(x-b) + c


azimuthData = np.deg2rad(azimuth)
coeffs, matcov = sc.curve_fit(cos_Wave, azimuthData, velocity, p0=[1, 0, 0])

plt.scatter(azimuthData, velocity)
nx = np.linspace(0, 2 * np.pi, 100)
plt.plot(nx, cos_Wave(nx, *coeffs))
plt.savefig("plot.png")
print(coeffs)

[ 6.63878549 1.03148322 -0.27674095]

enter image description here

David Jaimes
  • 331
  • 3
  • 7
  • Ah, I feel very dumb. I confused azimuth with elevation and I see where I went wrong now. It was a long Monday... Thanks for the help! – bgoudeau Jun 27 '17 at 16:54
  • 1
    @bgoudeau Are the values in `azimuth` supposed to be the 8 compass points (i.e. [0, 45, 90, 135, 180, 225, 270, 315])? If so, then the array `azimuth` should be created with, for example, `azimuth = np.linspace(0, 360, 8, endpoint=False)`. Currently in this answer, `azimuth` contains 0 *and* 360, and the spacing between the angles is 51.43. – Warren Weckesser Jun 29 '17 at 15:47
  • Yes, they are. That was a correction I made in my code after realizing what I was doing wrong after seeing the solution put forth by David. Thanks for the help Warren, sorry about not posting the full code as well. I am fairly sure that this one is solved. – bgoudeau Jun 29 '17 at 16:20