I have many individual datasets of $(x,y,z)$ coordinates. I'd like to see if the coordinates form a plane (and if they do, I'd then like to compute a unit normal to the plane). My original method was to simply fit to $z=ax+by+c$ via least-squares minimization. However, some of my data is nearly vertical in $z$ such that this breaks down. What is the best way to address this without specifically checking every dataset to see if the data forms a vertical plane?
Here is a toy dataset: coords = [3.64811 7.61531 9.05108; 3.53604 4.82801 9.05108; 3.53604 4.82801 6.34192; 3.64811 7.61531 6.34192]