I don't like the idea of someone believing me because I sound trustworthy, so let me waffle on a bit about the differences between two consecutive single-variable polynomial fits, versus a single bi-variate polynomial fit. For simplicity, I'll limit to cubic (third degree) polynomials, only because it keeps the number of coefficients down.
Let's say \$x\$ is the pressure reading, and \$y\$ is the temperature reading, directly from the hardware (ADC).
RTDs are almost linear, but not quite. Their resistance as a function of temperature is given by
Callendar–Van Dusen equation,
$$R(T) = R(0) \biggr( 1 + A T + B T^2 + (T - 100) C T^3 \biggr) = R_0 + R_1 T + R_2 T^2 + R_3 T^3 + R_4 T_4$$
where \$R_0 = R(0)\$, \$R_1 = R(0) A\$, \$R_2 = R(0) B\$, \$R_3 = -100 R(0) C\$, and \$R_4 = R(0) C\$; and the three coefficients, \$A\$, \$B\$, and \$C\$, are dependent on the actual sensor (material).
The above equation can be solved algebraically for \$T\$ when \$R(T)\$ is known, but the expression is absolutely horrible: nobody uses it. The three coefficients are very close to zero, so \$R(T)\$ is almost linear, \$R(T) \approx R_C (T - T_0)\$. Therefore, \$T(R)\$ is almost linear too, and we can approximate it using a polynomial also. Typically, a quartic (fourth degree) polynomial is used.
To measure the resistance, we use an
analog to digital converter, ADC, which converts either voltage or current to a reading. Typically – for example in the ubiquitous MAX31865 RTD-to-digital converter –, the RTD is fed from a precise current source, and the voltage (difference, in
four-terminal configuration) converted to a reading using a linear ADC.
To measure the
pressure, there are a number of completely different types of sensors. Piezoresistive strain gauges and metal strain gauges have a pretty linear response in their stable measurement region, so we can assume the pressure reading is more or less linear function of the actual pressure. Depending on the actual sensor, some compensation due to nonlinearity may be needed, so often a polynomial fit (with zeroth and first degree terms being "large", and the higher degree terms very small; so "almost linear") is used.
Every one of these is an approximation, fitted to the measured parameters. We don't have any specific model, except for, uh,
"this is almost linear". (We do know that for RTDs, second and higher-degree terms are all negative for most common materials used in RTDs, for physical reasons actually.) So, the sanity check for the models is to examine a few (say, two to five) different temperatures, and see how it affects a specific (constant) pressure; and a few different pressures, to see how temperature variation affects the pressure reading. If they look to be within expected values, the overall approximation is valid.
So, what's the difference between applying consecutive polynomial fits, and a single multi-variate fit? Cross terms.
Consider the situation when you have two different readings, \$x\$ and \$y\$, and you have a function \$f(x, y)\$ to evaluate, but both \$x\$ and \$y\$ are only almost linear. You can calculate \$f\bigr(X(x), Y(y)\bigr)\$ where \$X(x)\$ is the compensated \$x\$ and \$Y(y)\$ is the compensated \$y\$. This is perfectly okay to do, if you have a reason and reliable fits for \$X(x)\$ and \$Y(y)\$.
In this case, \$f(x, y) = x + g(y)\$, i.e. temperature \$y\$ is only used to give a correction \$g(y)\$ to the pressure \$x\$, to estimate the actual pressure.
The pressure compensation is only dependent on temperature, not pressure. In fact, we could write
$$f(x, y) = x + G_1 (y - y_0) + G_2 (y - y_0)^2 + G_3 (y - y_0)^3 + G_4 (y - y_0)^4$$
where \$y_0\$ is the temperature at which the pressure readings are correct with zero compensation. See how there are no terms multiplying both \$x\$ and \$y\$? It means this model cannot vary the pressure compensation based on pressure
and temperature, only on temperature. I do not believe this is physically sound.
Consider the case where the estimated actual pressure \$f(x, y)\$ is a bivariate (temperature reading \$y\$, pressure reading \$x\$) polynomial of degree four:
$$\begin{aligned}
f(x,y) &= F_{0 0} + F_{0 1} y + F_{0 2} y^2 + F_{0 3} y^3 + F_{0 4} y^4 + \\
~ & ~ ~ F_{1 0} x + F_{1 1} x y + F_{1 2} x y^2 + F_{1 3} x y^3 + F_{1 4} x y^4 + \\
~ & ~ ~ F_{2 0} x^2 + F_{2 1} x^2 y + F_{2 2} x^2 y^2 + F_{2 3} x^2 y^3 + F_{2 4} x^2 y^4 + \\
~ & ~ ~ F_{3 0} x^3 + F_{3 1} x^3 y + F_{3 2} x^3 y^2 + F_{3 3} x^3 y^3 + F_{3 4} x^3 y^4 + \\
~ & ~ ~ F_{4 0} x^4 + F_{4 1} x^4 y + F_{4 2} x^4 y^2 + F_{4 3} x^4 y^3 + F_{4 4} x^4 y^4 \\
\end{aligned}$$
Many of the cross terms (\$F_{i j}\$ where \$i \ne 0\$, \$j \ne 0\$) are zero or very close to zero. They describe the change in the pressure compensation as a function of
both pressure and temperature.
I believe these cross terms are needed, because we do not know exactly the pressure sensor behaviour as a function of temperature at different pressures (i.e, \$f(x, y)\$ in exact form!), but there is reason to believe that while it is likely small, there is some, due to mechanical reasons alone. (For example, thermal expansion of the strain gauge, and the fixture where the strain gauge is located, if the pressure sensor uses the strain gauge.)
If mechanical, these dependencies are complex, but small. We have about zero hope of deriving a numerical description from first principles – we
can compute it, it is just way too complicated to be worth the effort.
Instead, we do a set of measurements at known temperatures and pressures. If we draw a diagram with pressure on one axis, and temperature on the other axis, using a polynomial fit to a number of samples gives us an interpolated estimate within the convex hull of the samples, and handwavy extrapolated estimate outside the convex hull. (Outside the convex hull we are basically extrapolating the system behaviour, that's the hand-waving bit.) The denser the samples, the more reliable the estimate in that region.
If we only do measurements along each axis (i.e., one set of samples roughly horizontally in one line, and another set of samples roughly vertically on one line), we have no information for the cross terms, and fitting the measurements separately in their single-variate polynomials makes more sense.
(I like to use Gnuplot for fitting, because it allows/requires one to provide the function used for the fit, and the variables and coefficients used in the fitting. I haven't used Excel/LibreOffice Calc for this, but that does not mean it's not a valid tool also; I'm just unfamiliar with that one. I do know others use it, although Matlab/Octave is even more popular among scientists doing numerical work.)