I find it a bit easier to work with homogeneous coordinates and matrices for this. In the following, I use lower-case bold letters to represent homogenous vectors and a tilde to indicate their corresponding inhomogeneous coordinate tuples.
Your general conic equation can be written as $$\mathbf x^TC\mathbf x = \begin{bmatrix}x&y&1\end{bmatrix} \begin{bmatrix} a&\frac b2&\frac d2 \\ \frac b2 & c & \frac e2 \\ \frac d2&\frac e2&f \end{bmatrix} \begin{bmatrix}x\\y\\1\end{bmatrix} = 0.$$ If you have an invertible point transformation $\mathbf x' = M\mathbf x$, then $$\mathbf x^TC\mathbf x = (M^{-1}\mathbf x')^TC(M^{-1}\mathbf x') = \mathbf x'^T(M^{-T}CM^{-1})\mathbf x',$$ so the conic’s matrix transforms as $C' = M^{-T}CM^{-1}$.
The invertible affine transformation $\tilde{\mathbf x} = A\tilde{\mathbf x}'+\tilde{\mathbf t}$ (note that I’m inverting your transformation for simplicity) can be represented by the $3\times 3$ matrix $$M^{-1} = \left[\begin{array}{c|c} A & \tilde{\mathbf t} \\ \hline \mathbf 0^T & 1\end{array}\right]$$ so that $\mathbf x = M^{-1}\mathbf x'$. Writing $C$ in block form as $$C = \left[\begin{array}{c|c} Q & \tilde{\mathbf b} \\ \hline \tilde{\mathbf b}^T & f \end{array}\right],$$ which corresponds to the form $\tilde{\mathbf x}^TQ\tilde{\mathbf x}+2\tilde{\mathbf b}^T\tilde{\mathbf x}+f = 0$ of the general conic equation, we then have $$M^{-T}CM^{-1} = \left[\begin{array}{c|c} A^T & \mathbf 0 \\ \hline \tilde{\mathbf t}^T & 1\end{array}\right] \left[\begin{array}{c|c} Q & \tilde{\mathbf b} \\ \hline \tilde{\mathbf b}^T & f \end{array}\right] \left[\begin{array}{c|c} A & \tilde{\mathbf t} \\ \hline \mathbf 0^T & 1\end{array}\right]
= \left[\begin{array}{c|c} A^TQA & A^T(Q\tilde{\mathbf t}+\tilde{\mathbf b}) \\ \hline (Q\tilde{\mathbf t}+\tilde{\mathbf b})^T A & \mathbf t^TC\mathbf t \end{array}\right].$$ There are several things to note here: the quadratic part of the equation represented by $Q$ is only affected by the linear part of the affine transformation; the new constant term is the left-hand side of the untransformed equation evaluated at $\mathbf t$, the translation part of the transformation; and the signs of $\det C$ and $\det Q$ are preserved by affine transformations. If you examine the signs of these determinants for a canonical example of each the types of conics that you’ve listed, you’ll find that the combination of signs determines the type of conic—opposite signs for an ellipse, $\det Q=0$ and $\det C\ne0$ for a parabola, $\det C=0$ and $\det Q\lt0$ for a pair of intersecting lines, and so on. In fact, $\det Q$ is a multiple of the discriminant of the conic equation. It’s a bit tedious, but not very difficult, to work out appropriate transformations to bring an arbitrary conic into one of these canonical forms. (You might need to use the fact that a real symmetric matrix is orthogonally diagonalizable to prove existence of $A$ in all cases.) For example, $\det Q\ne0$ for any central conic, so the linear terms in the equation of a central conic can be eliminated by choosing $\tilde{\mathbf t} = -Q^{-1}\tilde{\mathbf b}$.
[Transferring the discussion in comments to an answer.]
Under a general affine transformation, the image of the center of an ellipse is the center of the ellipse’s image, so that’s easily computed. Unfortunately, the principal axes aren’t usually mapped to the axes of the image, although their images are conjugate diameters of the resulting ellipse. There are a couple of approaches that you might take to determine the semiaxis lengths and rotation angle of the image.
Writing the general conic equation $Ax^2+Bxy+Cy^2+Dx+Ey+F=0$ in matrix form $$\mathbf x^TQ\mathbf x = \begin{bmatrix}x&y&1\end{bmatrix} \begin{bmatrix}A&\frac B2&\frac D2\\\frac B2&C&\frac E2\\\frac D2&\frac E2&F\end{bmatrix} \begin{bmatrix}x\\y\\1\end{bmatrix} =0, \tag1$$ under the point transformation $\mathbf x' = M\mathbf x$ the matrix $Q$ transforms into $M^{-T}QM^{-1}$. (You can easily derive this by substituting for $\mathbf x$ in (1) and rearranging.) So, you could compute this transformation, extract the new coefficients of the conic equation from the result, and then use standard methods and formulas to extract the parameters that you’re interested in.
You can save yourself some work by looking only at the quadratic part of the equation: the shape of the ellipse is only affected by the linear part of the transformation, so you need only focus on the ellipse with equation $${(x\cos\theta+y\sin\theta)^2\over a^2}+{(x\sin\theta-y\cos\theta)^2\over b^2}=1, \tag 2$$ i.e., decompose the affine transformation into the composition of a translation by $(-h,-k)$, a linear transformation, and then another translation. (Note that I’ve corrected the sign error in your original equation.) In matrix form, we have $(x,y)Q(x,y)^T=1$ and the ellipse’s matrix again transforms as $M^{-T}QM^{-1}$, except that now we’re working with $2\times2$ matrices. Apply the transformation and extract the rotation and semiaxis lengths from the resulting coefficients.
One could, with a lot of tedious work or assistance of a symbolic algebra program, work out a general formula for the resulting parameters, but it’s going to be very ugly. If you’re planning to code this, it’s probably best to take it step-by-step instead of trying to transcribe some complex formula.
There are a couple of alternative approaches that might be interesting, though they don’t lead to simple general formulas, either, and might not be any less work than the above. First notice that you can parameterize your original ellipse as $\mathbf c + \mathbf u\cos t + \mathbf v\sin t$, where $\mathbf c = (h,k)$, $\mathbf u = a(\cos\theta,\sin\theta)$ and $\mathbf v=b(-\sin\theta,\cos\theta)$. Focusing again on the linear part of the transformation, the new ellipse is a translate of $\mathbf r(t) = M\mathbf u\cos t+M\mathbf v\sin t$. The tangent to an ellipse at a vertex is orthogonal to the corresponding principal axis, so a way to find these vertices is to solve $\mathbf r(t)\cdot\mathbf r'(t)=0$ for $t$. Unless you ended up with a circle, there will be two pairs of diametrically opposite solutions. The semiaxis lengths are then just the two lengths $\|\mathbf r(t)\|$ and the rotation angle is also easily computed.
Another possible approach takes advantage of the fact that the original ellipse (again after translation back to the origin) is itself obtained from the unit circle by scaling and rotation, so that this ellipse when further transformed is obtained from the unit circle via some composite transformation given by the matrix $M'=MRS$, with $$R = \begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix} \text{ and } S = \begin{bmatrix}a&0\\0&b\end{bmatrix}.$$ The SVD of $M'$ refactors it into the product of a rotation, a diagonal matrix (i.e., scaling along orthogonal directions) and another rotation. The half-axis lengths of the ellipse are therefore given by the central diagonal matrix, and, since the first rotation leaves the unit circle unchanged, the ellipse’s principal axis directions are the columns of the left-hand rotation matrix in the SVD. The rotation angle is easily recovered from them.
Best Answer
To keep it understandable albeit inelegant I'll do the passage in two steps
First the translation
$x^2 - xy + y^2 - 3y -1 = 0$
we look for a new centre $(h;\;k)$ so we substitute $x=x'+h;\;y=y'+k$
$-(h+x) (k+y)+(h+x)^2+(k+y)^2-3 (k+y)-1=0$
$x'^2-x' y'+y^2+x' (2 h-k)+y' (-h+2 k-3) +h^2-h k+k^2-3 k-1=0$
To have no first degree terms we put $2h-k=0;\;-h+2 k-3=0$ which gives
$h=1;\;k=2$ and we plug these values in the previous equation
$x'^2 - x' y' + y'^2=4$
Now we want to get rid of the $x'y'$ term. To do so we have to rotate the axis using these equations
$\begin{aligned}x'&=X\cos \theta -Y\sin \theta \\y'&=X\sin \theta +Y\cos \theta \end{aligned}$
$(X \sin \theta+Y \cos \theta)^2-(X \cos \theta-Y \sin \theta) (X \sin \theta+Y \cos \theta)+(X \cos \theta-Y \sin \theta)^2=4$
collecting terms
$X^2 \left(\sin ^2\theta+\cos ^2\theta-\sin \theta \cos \theta\right)+X Y \left(\sin ^2\theta-\cos ^2\theta\right)+Y^2 \left(\sin ^2\theta+\cos ^2\theta+\sin \theta \cos \theta\right)=4$
as we want the term $XY$ off we set $\sin ^2\theta-\cos ^2\theta=0$
$\tan^2\theta=1\to \theta=\pm \dfrac{\pi}{4}$
if we want the major axis of the ellipse to be horizontal we choose $\theta=\dfrac{\pi}{4}$ and substitute this value in the last equation
$\dfrac{X^2}{8}+\dfrac{3 Y^2}{8}=1$
the equations of the roto-translation altogether are
$\begin{aligned}x&=X\cos \frac{\pi}{4}-Y\sin \frac{\pi}{4} +1\\y&=X\sin \frac{\pi}{4}+Y\cos \frac{\pi}{4} +2\end{aligned}$
or
$\begin{aligned}x&=\frac{\sqrt 2}{2}(X-Y) +1\\y&=\frac{\sqrt 2}{2}(X+Y) +2\end{aligned}$
hope this helps