As I went around the little jogging track, I realized that I wanted to know just how much of the hill would be covered by Kistahöjden[a], the new residential area that is to be built there.

2012-12-03 12:10

In order to do that, I had to fit the map from Kistahöjden's web site[b] onto Google Earth[c]. This turned out to be quite difficult. I pulled and stretched and rotated the image, but couldn't quite get it to line up. Finally I figured, why not just whip up some code that will compute the transform parameters for me? This is the result.

​1. The Application

Here it is, for the impatient:
 Image Width px Height px Point 1 Latitude deg Longitude deg X px Y px Point 2 Latitude deg Longitude deg X px Y px Rotation deg Center Latitude deg Longitude deg Unrotated Bounds Suitable for use in Google Earth along with the rotation value North deg East deg South deg West deg Rotated Corners Bottom-Left deg deg Top-Right deg deg

​2. Instructions

The theory behind the application is that given two pairs of reference points, with one of the pair in the overlay image and one in the map, there is only one unique way to translate and rotate the overlay image so that the overlay points in the pairs coincide with the corresponding map points.

 ! The image we will try to overlay is on this page[d]. I've tried to contact the intellectual property holders of the map, but to no avail. Since the Swedish intellectual property laws make no exception for "quoting" maps, I can't provide any illustrative images here. You should, however, be able to download the map yourself and use it in Google Earth.

​2.1. Choosing Reference Points

We start by choosing the reference points. Good points are for example centers of crossings or roundabouts. Then we measure the image and map coordinates of the reference points.

​2.2. Using the Values

With these values, we just have to add the width and height of the overlay image and we're done. To save you the hassle of copying and pasting the values:

​3. Result

The algorithm can be used to create an image overlay on a Google Maps map:

``````var overlay = monochrome.maps.ImageOverlay.fromReferencePoints (
{
lat: 59.412744,
lon: 17.932549,
x: 270,
y: 582
},

{
lat: 59.411911,
lon: 17.941253,
x: 845,
y: 190
},

{
width: 960,
height: 683
}
);``````

In the overlay below, I've replaced the actual map with two big red spots where the roundabouts were due to IP issues. You should, however, be able to see that the red spots match the roundabouts.

​4. The Math Behind

The basic idea is to decompose the transformation from image coordinates to world ("map" or latitude / longitude) coordinates into three steps:

1. A rotation that aligns the line between the two reference points in the image with the line between the two reference points in the map.

2. A translation that moves the image so that the first reference point pair coincides.

3. A scaling that aligns the second reference point pair.

I will refer to the two reference point pairs as  p_(1)  and  p_(2) . Each point pair has an x and a y coordinate - for the point in the image, as well as a latitude and a longitude coordinate - for the point on the map. These four numbers will be called  p_(n,x) ,  p_(n,y) ,  p_(n,lat) ,  p_(n,lon)  for  n = 1,2 .

First, let's compute the deltas:

(eq.1)

Delta_(x) = p_(2,x)-p_(1,x)

Delta_(y) = p_(2,y)-p_(1,y)

Delta_(lat) = p_(2,lat)-p_(1,lat)

Delta_(lon) = p_(2,lon)-p_(1,lon)

​4.1. Rotation

We define the "width" of a degree of latitude or longitude as the distance along the surface between two parallells (lines of constant latitude) or meridians (lines of constant longitude). As we move north or south from the equator, the width of each longitude degree shrinks, and reaches zero at the poles. When computing angles in lat/long space, we need to compensate for that. Fortunately, we can compute an "average latitude per longitude ratio" for the overlay and not be far off. For a latitude  lat , the ratio of the surface length of a latitude degree to longitude degree is:

(eq.2)

lat:lon = 1:cos(lat)

That is, at the equator, one degree of latitude is as wide as one degree of longitude, but halfway to the pole, one degree of longitude is only ~70% of the width of one degree of latitude.

A suitable number to use for  lat  is simply the average latitude of the reference points:  (p_(1,lat) + p_(2,lat)) / 2 . Once we have this correction factor we can assume that lat/long space is a cartesian grid. Finding the rotation is then easy: The angle of the line between two points relative to the positive x-axis is  tan^-1(Delta_(y),Delta_(x)) , so we just subtract the resulting angle from one pair of points from the other:

(eq.3)

theta_(c) = cos((p_(1,lat)+p_(2,lat))/2)

alpha_(map) = tan^-1(Delta_(lat), Delta_(lon)theta_(c))

alpha_(img) = tan^-1(Delta_(y),Delta_(x))

r = alpha_(img) - alpha_(map)

​4.2. Translation

Finding the translation is also easy: Since we only have to satisfy the constraint that  p_(1,x),p_(1,y)  shall correspond to  p_(1,lat),p_(1,lon) , a transformation is:

(eq.4)

f(x,y) -> [ lon', lat' ]

lon' = x-p_(1,x)+p(1,lon)

lat' = y-p_(1,y)+p(1,lat)

​4.3. Scale

Finally, we need to find the scaling factors. One for number of latitude degrees per pixel, and one for the number of longitude degrees per pixel:  s_(lat)  and  s_(lon) . This turned out to be the hardest part, because we have four cases - three, if we ignore the degenerate case. It is obvious that if  Delta_(x)  and  Delta_(y)  are zero, which implies that  Delta_(lat)  and  Delta_(lon)  also are zero, then no meaningful scale value can be computed, as we then only have a single reference point. However, it is perfectly possible for one of each pair to be zero. This gives us two more cases. We will, however, start with the last case: when all four deltas are non-zero.

Imagine standing at  p_(1) , on an image that is rotated and translated. You are looking down the positive x-axis of the image. Then you walk  Delta_(x)  pixels forward. As you do so, you move  s_(lon)cos(r)Delta_(x)  degrees of longitude and  s_(lat)sin(r)Delta_(x)  degrees of latitude. Then you turn 90 degrees and walk along the y-axis of the image. You then move  s_(lon)cos(r + pi/2)Delta_(y)  degrees of longitude and  s_(lat)sin(r + pi/2)Delta_(y)  degrees of latitude.

In total, then, you have moved  s_(lon)cos(r)Delta_(x) + s_(lon)cos(r + pi/2)Delta_(y)  degrees of longitude and  s_(lat)sin(r)Delta_(x) + s_(lat)sin(r + pi/2)Delta_(y)  degrees of latitude. But if our values for  s_(lat)  and  s_(lon)  are correct, these distances along the latitude and longitude axes should be equal to  Delta_(lat)  and  Delta_(lon) ! We therefore set:

(eq.5)

r_(2) = r+pi/2

Delta_(lon) = s_(lon)cos(r)Delta_(x) + s_(lon)cos(r_(2))Delta_(y)

Delta_(lat) = s_(lat)sin(r)Delta_(x) + s_(lat)sin(r_(2))Delta_(y)

...and solve for  s_(lat)  and  s_(lon) :

(eq.6)

s_(lon) = Delta_(lon) / (cos(r)Delta_(x) + cos(r_(2))Delta_(y))

s_(lat) = Delta_(lat) / (sin(r)Delta_(x) + sin(r_(2))Delta_(y))

​4.3.1. The Zero-delta Cases

For the cases when  Delta_(lat)  or  Delta_(lon)  are zero, we have to rely on the average ratio of latitude to longitude width. For these two cases, then, we compute the scale factor for the non-zero delta, and multiply or divide as needed with the average latitude per longitude ratio to get the other scale factor.

​4.4. The Transformation

Putting it all together, we find that the transformation is:

(eq.7)

f([x,y]) -> [lat' , lon']

delta_(x) = x-p_(1,x)

delta_(y) = y-p_(1,y)

lat' = p_(1,lat) + sin(r)*delta_(x)*s_(lat) + sin(r_(2))*delta_(y)*s_(lat)

lon' = p_(1,lon) + cos(r)*delta_(x)*s_(lon) + cos(r_(2))*delta_(y)*s_(lon)

We can now compute the latitude and longitude of any point in the image. Placing the image in the map is now trivial: We compute the latitude and longitude for the lower left corner, use the `fromLatLngToDivPixel` method to figure out screen-space coordinates, apply CSS width and height values to scale and a CSS transformation to rotate around the lower left corner.

​5. Source and Docs

2012-12-22, updated 2013-03-09