davenquinn / Attitude

Attitude is for fitting the orientation of planes!
MIT License
11 stars 2 forks source link

PCAOrientation covariance_matrix incorrect for plotting vs ReconstructedPlane #10

Open AndrewAnnex opened 2 years ago

AndrewAnnex commented 2 years ago

While attempting to plot some real data with the stereonet or polar plotting methods I found that the error ellipses plotted did not make sense relative to the fitted data results. I believe I isolated the issue to the default covariance_matrix returned by PCAOrientation relative to the one returned by converting the same fit to a ReconstructedPlane object.

for example a PCAOrientation fitted as:

Orientation:: strike: 187.28 dip: 14.42
      error:: min: 1.99 max: 24.30 rake: 92.55

plots as the following using the object directly using a polar plot

image

constructing a ReconstructedPlane and modifying the following attributes using:

rp = attitude.ReconstructedPlane(*pcafit.strike_dip_rake(),*pcafit.angular_errors()[::-1])
def _sd(): 
    return rp.strike_dip_rake()[0:2]
rp.strike_dip = _sd #necessary for plotting
rp.azimuth = rp.dip_direction()[0] #necessary for plotting

results in the following plot using the same code

image

The only relevant difference I could find is the returned covariance matrix which for the pcafit object is:

array([[1.68404678e+07, 0.00000000e+00, 0.00000000e+00],
        [0.00000000e+00, 2.65137967e+05, 0.00000000e+00],
        [0.00000000e+00, 0.00000000e+00, 1.10673938e+02]]),

and for the rp object

 array([[94.01198358,  0.        ,  0.        ],
        [ 0.        , 23.41390221,  0.        ],
        [ 0.        ,  0.        ,  1.        ]]))

It's not clear why these are not equivalent

davenquinn commented 2 years ago

Yeah — the ReconstructedPlane has always had a few bugs. The math is not nearly as well-controlled as the primary forward fitting algorithms. It appears that the scaling on both error ellipses has some problems in this view, alas

These polar plotting functions in Python have always been a bit terrible. I'll try to give this another look later this week.