[I am working on researching a complete answer, this is as far as I've got so far, I will edit as I work out more of the solution.]
Further research has turned up this paper which provides at least part of the answer.
Briefly, ordinary least squares (OLS) does not consider uncertainties on either axis, and weighted least squares (WLS) only takes into account uncertainties in the predictor variable / dependent variable / on the y-axis. The above publication advocates the use of bivariate least squares (BLS) over orthogonal distance regression (ODR), but goes on to derive formulae for prediction intervals which appear to be compatible with ODR.
In this approach, each 'data point' $\left(x_{i}, y_{i}\right)$ is considered to be the result of multiple experiments or observations under ostensibly the same conditions, thus allowing the calculation of the variance for both the predictor and response variable pair.
The predictor ($x_{i}$) and response ($y_{i}$) variables are related by the following equation
$y_{i}=b_{0}+b_{1} x_{i}+e_{i}$
where $b_{0}$ and $b_{1}$ are the estimates of the intercept and slope of the true linear model, and $e_{i}$ is the $i$th residual error. The variance of $e_{i}$ is referred to as the weighting factor and is denoted as $w_{i}$ or $s_{e_{i}}^{2}$:
$w_{i}=s_{e_{i}}^{2}=s_{y_{i}}^{2}+b_{1}^{2} s_{x_{i}}^{2}-2 b_{1} \operatorname{cov}\left(x_{i}, y_{i}\right)$
where $s_{x_{i}}^{2}$ and $s_{y_{i}}^{2}$ are the experimental variances of point $i$, and $\operatorname{cov}\left(x_{i}, y_{i}\right)$ is the covariance between the measurements for each $\left(x_{i}, y_{i}\right)$ data pair.
The variance of the response variable $y_{0}$, being the mean of the $q$ observations performed at $x_{0}$, is given by:
$s_{y_{0}}^{2}=\left(\frac{1}{q}+ X _{0}^{ T }\left( X ^{ T } W ^{-1} X \right)^{-1} X _{0}+s_{x_{0}}^{2} b_{1}^{2}\right) s^{2}$
where $s^{2}=\frac{\sum_{i=1}^{n}\frac{\left(y_{i}-\hat{y}_{i}\right)^{2}}{w_{i}}}{n-2}$ is the estimate of the true experimental variance, $X_{0}$ is a two element column vector $\left|\begin{array}{l}
1 \\
x_{0}
\end{array}\right|$, and $X$ is an $n \times 2$ matrix for which the first column is a column of ones and the second column is formed by the $n$ values of $x$ corresponding to the experimental points; $W$ is an $n \times n$ diagonal matrix whose $i$th diagonal element is the weighting factor $w_{i}$, defined above.
The variance of the predictor mean value at a given observation $y_{0}$ is
$s_{x_{0}}^{2}=\left( Y _{0}^{ T }\left( Y ^{ T } W ^{\prime -1} Y \right)^{-1} Y _{0}+s_{y_{0}}^{2} \frac{1}{b_{1}^{2}}\right) s^{\prime 2}$
where $Y _{0}$ is $\left|\begin{array}{l}
1 \\
y_{0}
\end{array}\right|$, $Y$ is an $n \times 2$ matrix for which the first column is a column of ones and the second column is the values of $y$, $W^{\prime}$ is an $n \times n$ diagonal matrix whose $i$th diagonal element is the weighting factor $w_{i}^{\prime}=s_{x_{i}}^{2}+\frac{1}{b_{1}^{2}} s_{y_{i}}^{2}-2 \frac{1}{b_{1}} \operatorname{cov}\left(x_{i}, y_{i}\right)$, and $s^{\prime 2}$ the experimental error associated with the predictor variables, given by
$s^{\prime 2}=\frac{\sum_{i=1}^{n} \frac{\left(y_{i}-\hat{y}_{i}\right)^{2}}{w_{i}^{\prime}}}{n-2}$
The variance of the prediction of the predictor variable of a future sample at $y_{0}$, the mean of $q$ observations, is
$s_{x_{0}}^{2}=\left(\frac{1}{q}+ Y _{0}^{ T }\left( Y ^{ T } W ^{-1} Y \right)^{-1} Y _{0}+s_{y_{0}}^{2} \frac{1}{b_{1}^{2}}\right) s^{\prime 2}$
The prediction intervals of the response and predictor variables are then given by
$y_{0} \pm t_{\alpha, n-2} s_{y_{0}}$
$x_{0} \pm t_{\alpha, n-2} s_{x_{0}}$
where $t_{\alpha, n-2}$ is the $t$-value for the required significance level $\alpha$ and $n-2$ degrees of freedom.
TODO: Interface to scipy.odr output.
