Spherical to Cartesian coordinate ellipsoid overlap

Question

I have two geographic coordinates; latitude and longitude , separated by few meters. I need to draw an ellipsoid of same major and minor axes ,centered around the geographic co-ordinates

I used matplot lib and numpy to draw the ellipsoid. I am calculating the cartesian coordinates for the longitude using the formula given below To shift the ellipsoid to the lat and long co-rodinates, I am adding all the x,y,z array elements with catesian coordinates of longitude- x1, y1, z1

Is that the right way to shift the ellipsoid , so that if I draw another ellipsoid with a different geographic co-ordinate I will be able to visualize/calculate their overlap area

radius_of_world = 6371000
# Radii corresponding to the coefficients:
# http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node42.html ?
rx, ry, rz = 1/np.sqrt(coefs) # Is this correct for coefs gives as a,b,v

# Set of all spherical angles:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# Postion of the lattitude and logitude on surface of the earth// ellipsoidal earth
x1 = radius_of_world * cos(longitude) * cos(latitude)
y1 = radius_of_world * np.sin(longitude) * np.cos(latitude)
z1 = radius_of_world * np.sin(latitude)

print 'Point location  is', x1,y1,z1

# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
x = rx * np.outer(np.cos(u), np.sin(v)) 
y = ry * np.outer(np.sin(u), np.sin(v)) 
z = rz * np.outer(np.ones_like(u), np.cos(v))


# Plot:
#ax.plot_wireframe(x +x1 , y + y1, z + z1,  rstride=4, cstride=4, color=color)
ax.scatter(x +x1, y + y1, z + z1,color=color )

I am not sure if the calculation of rx, ry,rz is correct; I changed it slightly for seeing it in 2D and now I can see a and b coming with proper length

    Testing Elliptical Drawing in Spherical Corodinates

import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=(12,12))  # Square figure

def drawCoverage(coefs,latitude,longitude,color):


# Radii corresponding to the coefficients:
# http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node42.html ??
rx, ry = coefs # This is what I changed

# Set of all spherical angles:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# u = np.pi/2
# v = np.pi/2
# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
x = rx * np.outer(np.cos(u), np.sin(v)) 
y = ry * np.outer(np.sin(u), np.sin(v)) 

#ax.autoscale(enable=True)
# Plot:
#plt.axis([0,100,0,100])
plt.plot(x , y , 'bo' )


latitude,longitude= 13.04738626,77.61946793   
coefs = (373, 258)  # Coefficients in a0/c x**2 + a1/c y**2 + a2/c z**2 = 1 
drawCoverage(coefs,latitude,longitude,'b')

plt.show()

I changed the coef calculation to directly use a,b,c in the 3D model and now the axes are proper

# http://stackoverflow.com/questions/7819498/plotting-ellipsoid-with-matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=(12,12))  # Square figure
ax = fig.add_subplot(111, projection='3d')

def drawCoverage(coefs,latitude,longitude,color):

radius_of_world =  6372.8
# Radii corresponding to the coefficients:
# http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node42.html ??
rx, ry, rz = coefs

# Set of all spherical angles:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# Postion of the lattitude and logitude on surface of the earth// ellipsoidal earth
x1 = radius_of_world * cos(longitude) * cos(latitude)
y1 = radius_of_world * np.sin(longitude) * np.cos(latitude)
z1 = radius_of_world * np.sin(latitude)

print 'Point location  is', x1,y1,z1

# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
x = rx * np.outer(np.cos(u), np.sin(v)) 
y = ry * np.outer(np.sin(u), np.sin(v)) 
z = rz * np.outer(np.ones_like(u), np.cos(v))

#ax.autoscale(enable=True)
# Plot:
ax.plot_wireframe(x , y , z ,  rstride=4, cstride=4, color='b' )
ax.plot_wireframe(x + 111 , y + 111 , z + 111 ,  rstride=4, cstride=4, color='g')

#ax.scatter(x , y , z , color='b' )
#ax.scatter(x +111 , y +111 , z +111 , color='g')

#ax.set_xlim3d(0, 500)
#ax.set_ylim3d(0, 500)
#ax.set_zlim3d(0, 500)

# Adjustment of the axes, so that they all have the same span:
#max_radius = .25
#for axis in 'xyz':
#    getattr(ax, 'set_{}lim'.format(axis))((-max_radius, max_radius))



latitude,longitude= 13.04738626,77.61946793   
coefs = (373, 258, 258)  # Coefficients in a0/c x**2 + a1/c y**2 + a2/c z**2 = 1 
drawCoverage(coefs,latitude,longitude,'b')

latitude,longitude = 13.04638626, 77.61951605

# Distance between two lat.long is  0.1113 km
# http://www.movable-type.co.uk/scripts/latlong.html
#coefs = (258, 258, 373)
#drawCoverage(coefs,latitude,longitude,'g')
plt.show()

So what I did in the second code snippet is to calculate the distance between the two geographic locations using the Haversine formula ( which give distance as 111 meters approximately. I added this 111 to the second sphere for all x ,y and z points, so as to separate it out and see the overlap of the 'spheres of influence'. Is this approach correct

Answer 1

I guess this is the way

# drawing an ellipsoid - Oblate spheriod to show overlap for two geo points
at two altitudes

# http://stackoverflow.com/questions/7819498/plotting-ellipsoid-with-matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=(12,12))  # Square figure
ax = fig.add_subplot(111, projection='3d')

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

  R = 6372.8 # Earth radius in kilometers

  dLat = radians(lat2 - lat1)
  dLon = radians(lon2 - lon1)
  lat1 = radians(lat1)
  lat2 = radians(lat2)

  a = np.sin(dLat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dLon/2)**2
  c = 2*np.arcsin(np.sqrt(a))

  return R * c

def drawCoverage(coefs,horizontal_distance,vertical_distance):
    """Coverage Overalp of two cells sperated by horizotal distance between two geographic co-ordinates and
    vertical distance - height
    """

    # Radii corresponding to the coefficients:
    # http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node42.html ??
    rx, ry, rz = coefs

    # Set of all spherical angles:
    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)

    # Cartesian coordinates that correspond to the spherical angles:
    # (this is the equation of an ellipsoid):
    x = rx * np.outer(np.cos(u), np.sin(v)) 
    y = ry * np.outer(np.sin(u), np.sin(v)) 
    z = rz * np.outer(np.ones_like(u), np.cos(v))

    # Plot:
    ax.plot_wireframe(x, y, z ,  rstride=4, cstride=4, color='b' )
    ax.plot_wireframe(x +horizontal_distance, y, z + vertical_distance,rstride=4, cstride=4, color='g')

    latitude1,longitude1 = 13.04738626,77.61946793   
    latitude2,longitude2 = 13.04638626, 77.61951605
    dist = haversine(latitude2,longitude2,latitude1,longitude1)
    height_in_meters = 258*2                                 
    coefs = (373, 258, 258)  # 373 major axis, 258 miniir axis - oblate spheroid,roate along x
    drawCoverage(coefs,0,height_in_meters)# for visual validity using just height

plt.show()