Skew map to match aerial image
Well you could try to estimate a transformation matrix that maps between the coordinates in the two images. The idea is that

\begin{equation}
\begin{bmatrix}
x_a\\
y_a
\end{bmatrix}
=
\begin{bmatrix}
A_{xx} &amp; A_{yx}\\
A_{xy} &amp; A_{yy}
\end{bmatrix}
\begin{bmatrix}
x_m\\
y_m
\end{bmatrix}
\end{equation}

where $x_a, y_a$ are the coordinates (in pixels for example) of the aerial image and $x_m, y_m$ are the coordinates of the photograph. Writing this out leads to the system of equations:

\begin{equation}
x_a = A_{xx}x_m + A_{yx}y_m
\end{equation}
\begin{equation}
y_a = A_{xy}x_m + A_{yy}y_m
\end{equation}

That way you can represent scaling, rotation and shearing operations between the images. You need at least four points per image to solve that system but I recommend more and then calculate the least squares solution using linear regression. This will also give you an indication on how well everything works. If there are large residual errors after the fit, then the linear transformation might not be sufficient.

If you then want to overlay the images, just apply the transformation matrix to your map and the result should be an image that has the reference coordinates on top of each other and the areas in between are linearly transformed accordingly.

There are many things that can go wrong here. For once, it only works for roughly flat terrain  (even though that effect becomes weaker if the plane is further away). Also, if there is lens distortions in the photograph, that will affect the accuracy of the fit (usually close to the image borders). And finally, map inaccuracies might throw you off.

In the end, it's just a matter of trying. If the 2D-transformation model is not capable of delivering accurate results, you need to model more effects. But estimating those will probably also require more referene points.

One other idea would be to put the image on a plane (like a texture) and place the plane in a 3D-space and match it with the plane of the photograph. That would end you up with a 3-dimensional transformation matrix for this plane. Might work as well.