I’m not sure that there’s any particular coordinate system for $P$ that could be called “canonical.” If you have the defining equations for $P$, then you can build a coordinate system directly from them. For any $p,q\in P$, $p-q$ lies in the solution space $P'$ of the related homogeneous system of linear equations $\sum_iA_{1,i}x_k=0,\dots,\sum_iA_{k,i}x_i=0$, which is parallel to $P$ but passes through the origin. So, compute a basis for the kernel (null space) of the coefficient matrix $A$ of this system. If you need an orthonormal basis, apply the Gram-Schmidt process. If you only need it to be orthogonal, you can leave the resulting vectors unnormalized. Choosing any point $p_0\in P$ then gives you a coordinate system on $P$.
If, on the other hand, you want to build the coordinate system from some set of $k$ points, then the basic process you describe certainly works. Choose one of them, say $p_1$, to be the origin. The first basis vector is then $v_1=p_2-p_1$, normalized or not, as you like. The next basis vector will be the orthogonal rejection of $(p_3-p_1)$ from $v_1$, i.e., $v_2=(p_3-p_1)-{(p_3-p_1)\cdot v_1\over v_1\cdot v_1}v_1$, again, normalize or not as you see fit. The next basis vector $v_3$ will be the orthogonal rejection of $(p_4-p_1)$ from the span of these vectors, and so on. At each stage, you’ll be computing $v_i=(p_{i+1}-p_1)-\sum_{j=1}^{i-1}{(p_{i+1}-p_1)\cdot v_j\over v_j\cdot v_j}v_j$. Continue until you have $k$ basis vectors. If you ever get zero, the given set of points has a dependency and doesn’t actually define a $k$-plane. If you’re simply starting from a set of points and don’t really care about the dimension of their $m$-plane that they lie in, just throw that one out and keep going until you’ve processed all of the points.
You’ve now got a basis $\{v_i\}$ for $P'$. Recall that, by definition, the coordinates of a vector are the coefficients in its expression as a linear combination of basis vectors. So, to find the $P$-coordinates of a point $p\in P$, solve the equation $p-p_0=\sum_{i=1}^ka_iv_i$ for the unknown coefficients $a_i$. If there’s no solution, then $p\notin P$, of course. This can be done en masse by forming the matrix $B$ that has the $v_i$ as its columns, augmenting it with all of the points shifted by $p_0$, and row-reducing. The $P$-coordinates will be in the first $k$ rows. Alternatively, form the pseudo-inverse $(B^TB)^{-1}B^T$ and apply it to the shifted points, which can also be done in bulk. If the basis is orthogonal, the pseudo-inverse is particularly easy to compute: the $i$th row is $v_i/(v_i\cdot v_i)$. A disadvantage of this latter method is that it will happily assign $P$-coordinates to points that aren’t in $P$ (which might be an advantage if your data aren’t exact).
For example, let $P$ be defined by the equations $$\begin{align}x_1+x_3+2x_4&=12\\2x_1-x_2-x_4&=-4\end{align}$$ and points on this plane $p_1=(1,2,3,4)^T$, $p_2=(2,4,2,4)^T$ and $p_3=(4,9,2,3)^T$. Then $$A=\pmatrix{1&0&1&2\\2&-1&0&-1}$$ with kernel spanned by $(1,2,-1,0)^T$ and $(2,5,0,-1)^T$, which is easily found by row-reduction. We keep the first of these vectors for our kernel basis and take the orthogonal rejection of the second for the other basis vector: $$(2,5,0,-1)^T-{(2,5,0,-1)^T\cdot(1,2,-1,0)^T\over(1,2,-1,0)^T\cdot(1,2,-1,0)^T}(1,2,-1,0)^T=(0,1,2,-1)^T.$$ Taking $p_1$ as our origin, we form the matrix $$\left(\begin{array}{rr|rrr}1&0&1&3&2\\2&1&2&7&2\\-1&2&-1&-1&-1\\0&-1&0&-1&0\end{array}\right)$$ (I’ve thrown in a bogus point to illustrate error detection). Row-reduction yields $$\left(\begin{array}{rr|rrr}1&0&1&3&2\\0&1&0&1&-2\\0&0&0&0&5\\0&0&0&0&-2\end{array}\right)$$ and so the $P$-coordinates of $p_2$ are $B(p_2-p_1)=(1,0)^T$ and the $P$-coordinates of $p_3$ are $(p_3-p_1)=(3,1)^T$. The last column shows correctly that the third point was bogus.
For the pseudo-inverse, we have $$(B^TB)^{-1}B^T=\frac16\pmatrix{1&2&-1&0\\0&1&2&-1}$$ and applying it to the above shifted points produces $$\pmatrix{1&3&\frac76\\0&1&0}.$$ The first two columns agree with the previous result, but for the point not on our plane we’ve produced the coordinates of the point on the plane closest to it.
Starting from the given points instead, we have $v_1=p_2-p_1=(1,2,-1,0)^T$. Next, we compute $p_3-p_1=(3,7,-1,-1)^T$, and its orthogonal rejection from $v_1$ is $$(3,7,-1,-1)^T-{(3,7,-1,-1)^T\cdot(1,2,-1,0)^T\over(1,2,-1,0)^T\cdot(1,2,-1,0)^T}(1,2,-1,0)^T=(0,1,2,-1)^T.$$ This is, not coincidentally because of the way I derived the points for this example, the same as the above basis for $\ker A$, so the rest of the calculation proceeds as before.
