source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/distance_lon_lat.py @ 55

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import string
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid, cm
from netCDF4 import Dataset
#import arctic_map # function to regrid coast limits
#import cartesian_grid_test # function to convert grid from polar to cartesian
import scipy.special
import ffgrid2
import map_ffgrid
from matplotlib import colors



def distance_on_unit_sphere(lat1, lon1, lat2, lon2):

    # Convert latitude and longitude to 
    # spherical coordinates in radians.
    degrees_to_radians = math.pi/180.0
        
    # phi = 90 - latitude
    phi1 = (90.0 - lat1) * degrees_to_radians
    phi2 = (90.0 - lat2) * degrees_to_radians
        
    # theta = longitude
    theta1 = lon1 * degrees_to_radians
    theta2 = lon2 * degrees_to_radians
        
    # Compute spherical distance from spherical coordinates.
        
    # For two locations in spherical coordinates 
    # (1, theta, phi) and (1, theta, phi)
    # cosine( arc length ) = 
    #    sin phi sin phi' cos(theta-theta') + cos phi cos phi'
    # distance = rho * arc length
    
    cos = (math.sin(phi1) * math.sin(phi2) * math.cos(theta1 - theta2) + 
           math.cos(phi1) * math.cos(phi2))
    arc = math.acos( cos )

    # Remember to multiply arc by the radius of the earth 
    # in your favorite set of units to get length.
    return arc