10

For a section of a project I'm working on at University, I'm trying to recreate the Earth using Python and use that to plot specific locations on the surface and plot circles in various orientations, so they align with satellite data I have to give a representation of an airplane's location at a given time from the data set.

I started off by simply plotting a wireframe and plotted the points on the wireframe that I needed (all to scale the Earth and its geographical locale).

The issue I'm having is when I plot the points on a sphere-ish object with an image of the Earth overlayed on top, the points disappear when rotating the sphere past a certain point. So, the initial question is: How can I stop them from disappearing?

Second of all; I can't seem to find any way to plot circles centered on the sphere - for example, a circle around the equator and then manipulate that same idea to plot circles on the surface of the sphere to give an image like below:

Circles on the Earth

I know this is from google maps, but I'm curious if this can be done in Python (I assume so).

My current code is:

    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    from itertools import product, combinations
    import PIL
    
    #Plot the Earth
    f = plt.figure(1, figsize=(13,13))
    ax = f.add_subplot(111, projection='3d')
    u, v = np.mgrid[0:2*np.pi:30j, 0:np.pi:20j]
    x=6371*np.cos(u)*np.sin(v)
    y=6371*np.sin(u)*np.sin(v)
    z=6371*np.cos(v)
    ax.plot_wireframe(x, y, z, color="b")
    
    #GES ground station @ Perth & AES @ KLIA
    
    ax.scatter([-2368.8],[4881.1],[-3342.0],color="r",s=100)
    ax.scatter([-1293.0],[6238.3],[303.5],color="k",s=100)
    
    #Load earthmap with PIL
    bm = PIL.Image.open('earthmap.jpg')
    #It's big, so I'll rescale it, convert to array, and divide by 256 to get RGB values that matplotlib accept 
    bm = np.array(bm.resize([d/3 for d in bm.size]))/256.
    
    #d/1 is normal size, anything else is smaller - faster loading time on Uni HPC
    
    #Coordinates of the image - don't know if this is entirely accurate, but probably close
    lons = np.linspace(-180, 180, bm.shape[1]) * np.pi/180 
    lats = np.linspace(-90, 90, bm.shape[0])[::-1] * np.pi/180 
    
    #Repeat code specifying face colours 
    
    x = np.outer(6371*np.cos(lons), np.cos(lats)).T
    y = np.outer(6371*np.sin(lons), np.cos(lats)).T
    z = np.outer(6371*np.ones(np.size(lons)), np.sin(lats)).T
    ax.plot_surface(x, y, z, rstride=4, cstride=4, facecolors = bm)
    
    plt.show()

If there is any way I can get it so the points stop disappearing and plot even just a circle on the equator, that'd be great!

Thanks!

コリン
  • 1,005
  • 1
  • 5
  • 26
Nyancoil
  • 103
  • 1
  • 7
  • I guess you are aware of [`basemap`](http://matplotlib.org/basemap/users/examples.html) but just checking in case you didn't come across it before? Not sure about the rotational capabilities. However, doing an image search I did come across this (rather poor) [example](https://enumap.wordpress.com/2012/06/05/plot-base-map-by-python/) that appears to show a basemap projection inside an interactive matplotlib window. – roganjosh Feb 06 '17 at 13:24
  • @roganjosh no I hadn't heard of that, I'll look into it. Although, getting the picture of the earth isn't an issue really, it's more getting circles and the points disappearing when the sphere is rotated to be looking at them. – Nyancoil Feb 06 '17 at 13:30
  • The thing with basemap though is that it's designed to get contours etc. plotted on the surface of the Earth. Because it (potentially) has all your requirements in one package it might be easier to get them to play nicely rather than troubleshooting conflicts between 2+ different packages e.g. PIL. I'd be interested to know either way if you find it has the capabilities you need :) – roganjosh Feb 06 '17 at 13:37
  • 3
    Better example illustrating rotation with custom points plotted on the surface: https://waterprogramming.wordpress.com/2013/11/02/python-makes-the-world-go-round/ – roganjosh Feb 06 '17 at 14:02
  • 2
    It seems that basemap has the functionality of hiding the points behind the sphere already integrated. If you want to stay with the pure matplotlib approach, you may look at [this question here](http://stackoverflow.com/questions/41699494/how-to-obscure-a-line-behind-a-surface-plot-in-matplotlib). – ImportanceOfBeingErnest Feb 06 '17 at 22:16
  • not sure if it's helpful per se, but fyi give a look to OpenGL library, you might find useful insights. cheers – Lorenzo Bassetti Jan 04 '23 at 09:38

1 Answers1

0

I believe the first part of your question has mostly been resolved. In regard to the second part of your question, you can try the following in general,

%matplotlib ipympl
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


# Creating an arbitrary sphere
f = plt.figure(1, figsize=(10, 10))
ax = f.add_subplot(111, projection='3d')
u, v = np.mgrid[0:2*np.pi:100j, 0:np.pi:100j]
x=np.cos(u)*np.sin(v)
y=np.sin(u)*np.sin(v)
z=np.cos(v)
ax.plot_surface(x, y, z, color="b", alpha=0.3) #low alpha value allows for better visualization of scatter points



def draw_circle_on_sphere(p:float, a:float, v:float):
    '''
        Parametric equation determined by the radius and angular positions (both polar and azimuthal relative to the z-axis) of the circle on the spherical surface
        Parameters:
            p (float): polar angle
            a (float): azimuthal angle
            v (float): radius is controlled by sin(v)
            
        Returns:
            Circular scatter points on a spherical surface
    '''
    
    u = np.mgrid[0:2*np.pi:30j]
    
    x = np.sin(v)*np.cos(p)*np.cos(a)*np.cos(u) + np.cos(v)*np.sin(p)*np.cos(a) - np.sin(v)*np.sin(a)*np.sin(u)
    y = np.sin(v)*np.cos(p)*np.sin(a)*np.cos(u) + np.cos(v)*np.sin(p)*np.sin(a) + np.sin(v)*np.cos(a)*np.sin(u)
    z = -np.sin(v)*np.sin(p)*np.cos(u) + np.cos(v)*np.cos(p)

    return ax.scatter(x, y, z, color="r")


_ = draw_circle_on_sphere(3*np.pi/2, np.pi/4, np.pi/4) # example

Result

View result

Jason Yu
  • 101
  • 3