Python for Great Circle Projections – Part 2

Python for Great Circle Projections – Part 2

This post will be a bit mathematically involved. Bear with me.

Projections

Let’s jump right into the math. The formula for projecting any point in latitude, longtitude to its (x,y) equivalent on a rectangular plane is:

x = λcos(φ1)
y = φ

λ is the longitude;
φ is the latitude;
φ1 are the standard parallels (north and south of the equator) where the scale of the projection is true;
x is the horizontal position along the map;
y is the vertical position along the map.

This means that for every point on a great circle, we can find and plot the equivalent point on an equirectangular surface.

For our usage, we will be assuming the standard parallels are zero. Hence we can directly map longitude to x-axis and latitude to y-axis.

Curves

First of all, we need to convert a roll and pitch values into a latitude and longitude pair. Since for any two points on a sphere there exists a unique great circle, we can use (0,0) and one other value of lat, lng derived from the roll and pitch to get two points which we can use to draw the great circle.

def getCurve(roll, pitch, inres):

    #generates a lat lon pair from the roll and pitch. lon1 is always zero.
    lat1 = pitch
    lon1 = 0

The code above takes in roll and pitch in degrees, so we can take the lat1 to be the pitch directly (since latitude is the degrees above horizon). Lon1 will be 0. We now have the first coordinate pair (pitch, 0). Now we need one more.
Now we will calculate the second lat,lon point in the direction of the roll.

lat2 = math.asin( math.sin(lat1)*math.cos(0.01) + math.cos(lat1)*math.sin(0.01)*math.cos(roll) )
lon2 = lon1 + math.atan2(math.sin(roll)*math.sin(0.01)*math.cos(lat1),math.cos(0.01)-math.sin(lat2)*math.sin(lat2))

Next we find the distance between the points and see how many such distances can fit on a unit sphere. We do this so we can extend the line between the two points around the globe.

coor = []

# find the distance between the two points on a unit sphere
d = (math.acos(math.sin(lat1) * math.sin(lat2) + math.cos(lat1) * math.cos(lat2) * math.cos(lon1 - lon2)))

#find the number of such distances in the circumference of the unit sphere to complete round circle
rg = (2 * math.pi) / d + 5
f = 0

Next, we calculate the points on this great circle according to the resolution provided to the function (inres). This is a bit complicated, and the formulas are taken from here.

    #resolution of points (100 points for the curve)
    res = rg / inres
    while f < rg:
        #calculate points on great circle using formula here: http://williams.best.vwh.net/avform.htm and http://math.stackexchange.com/a/384719
        A = math.sin((1 - f) * d) / math.sin(d)
        B = math.sin(f * d) / math.sin(d)
        x = A * math.cos(lat1) * math.cos(lon1) + B * math.cos(lat2) * math.cos(lon2)
        y = A * math.cos(lat1) * math.sin(lon1) + B * math.cos(lat2) * math.sin(lon2)
        z = A * math.sin(lat1) + B * math.sin(lat2)

        lat_f = math.atan2(z, math.sqrt(x*x + y*y))
        lon_f = math.atan2(y,x)

        # print the points and paste on gpxvisualiser dot com to generate a google earth compatible kml file to see great circle on the globe
        # print math.degrees(lat_f),',',math.degrees(lon_f)

        #project the points on an equirectangular surface. x-axis is lons and y-axis is lats. Hence scale for both are -pi -> pi and -pi/2 -> pi/2 respectively.
        x = lon_f
        y = lat_f

        #add points to point arrays for pyplot to use later
        thisCoor = Coor(x,y)
        coor.append(thisCoor)

        f += res

    #plot the graph
    #plotme function given below
    coor.sort(key=lambda x: x.x, reverse=True )
    if __name__ == "__main__":
        plotme(coor)
    else:
        return coor

class Coor:
    def __init__(self, x,y):
        self.x = x
        self.y = y

Once this is processed, you can plot the coor array using the matplotlib library.

def plotme(coor):
    xax = [o.x for o in coor]
    yax = [o.y for o in coor]

    plt.plot(xax,yax)
    x1,x2,y1,y2 = plt.axis()

    # set the axis limits
    plt.axis((-math.pi,math.pi,-math.pi/2,math.pi/2))

    plt.show()

You will get curves as follows:

these are pyplot curves for different values of roll and pitch

the above code is available as a python module here.

You can use this to check for roll and pitch of the camera, but this is not very useful if one has to do the comparisons for each frame in the video manually. Let’s look at how to do this automatically using OpenCv in the next post.