El propósito de esta entrada es doble. Por un lado, construiremos una representación tridimensional interactiva del globo terrestre utilizando las librería gráfica Plotly para Python, y finalmente sobre este modelo gráfico representaremos la órbita de un satélite artificial, en este caso la ISS, para lo que utilizaremos la librería Skyfield y el método de los TLE descrito en mi anterior post.
El método seguido para generar una representación 3D interactiva del globo terrestre está fuertemente inspirado en la primera parte de este artículo de Ryota Kiuchi.
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
import numpy as np
from netCDF4 import Dataset
from skyfield.api import load, EarthSatellite
from astropy import units as u
from astropy.constants import R_earth
Descarga del mapa terrestre con datos del relieve¶
El mapa que utilizaremos como base para nuestra representación es el mapa topográfico del relieve de la superficie terrestre publicado por el NOAA con el nombre de ETOPO1 Global Relief Model. Concretamente descargaremos el fichero descrito como "grid-registered: netCDF". Se trata de un fichero en formato netCDF (Network Common Data Form), que es un formato de fichero adecuado para el intercambio de datos científicos en formato array. Estos ficheros son contenedores de un conjunto de variables similares a los arrays multidimensionales de Python Numpy.
El fichero ETOPO1 ETOPO1 Global Relief Model grid-registration contiene las elevaciones del terreno en el centro de cada elemento de una cuadrícula que cubre un rango de [-180º, 180º] en longitud (variable 'x-range') y de [-90º, 90º] en latitud (variable 'y_range'), con una separación de las líneas de la cuadrícula de 1 arcominuto, o lo que es igual, 1/60 = 0.01666667 grados en cada dimensión (variable 'spacing').
El número de líneas en la cuadrícula será 360 * 60 + 1 = 21601 en la dimensión x (longitud), y de 180 x 60 + 1 = 10801 en la dimensión y (latitud). Estos valores son los indicados en la variable 'dimension'. El elemento adicional que estamos sumando en cada dimensión se debe a que, en la versión "grid-registered" del mapa hay una línea de celdas duplicadas a lo largo de la línea de longitud -180/180 y de los polos (latitud -90/90). Esto permite una mayor exactitud en las elevaciones, aunque también se puede utilizar la versión cell-registered de este mismo mapa que carece de esa línea adicional.
Finalmente, la variable 'z' contendrá la elevación media en el centro de cada elemento de la cuadrícula. Estos datos, expresados en metros, estarán recogidos en forma de array unidimensional con un tamaño de 21601 x 10801 = 233312401.
Una vez descargado el fichero en una carpeta de nuestro ordenador, lo leeremos con la librería Python netCDF4.
datos = Dataset("Datos/ETOPO1_Ice_g_gdal.grd", "r")
Vemos que las variables son accesibles mediante un diccionario Python:
print(datos.variables)
{'x_range': <class 'netCDF4._netCDF4.Variable'>
float64 x_range(side)
units: Latitude
unlimited dimensions:
current shape = (2,)
filling on, default _FillValue of 9.969209968386869e+36 used, 'y_range': <class 'netCDF4._netCDF4.Variable'>
float64 y_range(side)
units: Latitude
unlimited dimensions:
current shape = (2,)
filling on, default _FillValue of 9.969209968386869e+36 used, 'z_range': <class 'netCDF4._netCDF4.Variable'>
float64 z_range(side)
units: meters
unlimited dimensions:
current shape = (2,)
filling on, default _FillValue of 9.969209968386869e+36 used, 'spacing': <class 'netCDF4._netCDF4.Variable'>
float64 spacing(side)
unlimited dimensions:
current shape = (2,)
filling on, default _FillValue of 9.969209968386869e+36 used, 'dimension': <class 'netCDF4._netCDF4.Variable'>
int32 dimension(side)
unlimited dimensions:
current shape = (2,)
filling on, default _FillValue of -2147483647 used, 'z': <class 'netCDF4._netCDF4.Variable'>
int32 z(xysize)
scale_factor: 1.0
add_offset: 0.0
_FillValue: -2147483648
node_offset: 0
unlimited dimensions:
current shape = (233312401,)
filling on}
Vamos a examinar cómo se puede acceder a cualquiera de las variables. Por ejemplo la variable elevación "z". Al estar las variables en un diccionario accedemos a ellas indexando con su nombre entre corchetes. La variable contiene algunos metadatos adicionales sobre la misma. En particular vemos que su shape es el propio de un array unidimensional.
datos.variables['z']
<class 'netCDF4._netCDF4.Variable'>
int32 z(xysize)
scale_factor: 1.0
add_offset: 0.0
_FillValue: -2147483648
node_offset: 0
unlimited dimensions:
current shape = (233312401,)
filling on
Y podremos acceder a su contenido mediante el procedimiento habitual:
z = datos.variables['z'][:]
z
masked_array(data=[-4228., -4228., -4228., ..., 2745., 2745., 2745.],
mask=False,
fill_value=1e+20)
Es decir, estamos ante un array de Numpy de tipo 'masked'. Como la máscara es "False" ello significa que todos los datos del array son válidos. No hay ninguno enmascarado.
Vamos también a comprobar que las dimensiones de la cuadrícula son las indicadas en la descripción anterior:
dim_lon, dim_lat = datos.variables['dimension'][:]
dim_lon, dim_lat
(21601, 10801)
Primera representación de los datos (quick-and-dirty)¶
Vamos a hacer una primera representación gráfica rápida y sin preocuparnos del formato, solo para comprobar la bondad de los datos. En primer lugar convertiremos el array z en un array bidimensional, y como no necesitamos una resolución tan alta en nuestra imagen, haremos una reducción del muestreo tomando una separación entre pixeles que corresponda a un grado de separación angular, es decir, dividiremos por 60 la resolución original de un arcominuto. De esta manera reduciremos la ocupación de memoria y las representaciones gráficas serán mucho más rápidas.
Z = np.reshape(z, (dim_lat,dim_lon))
Z = Z[::60, ::60]
Z.shape
(181, 361)
px.imshow(Z, color_continuous_scale='gray')
Conversión a coordenadas cartesianas sobre la superficie esférica¶
El siguiente paso va a ser la conversión de las coordenadas de longitud y latitud a coordenadas cartesianas sobre la superficie esférica, con el fin de representar un objeto Surface de Plotly. Para ello utilizaremos las fórmulas trigonométricas habituales, pero previamente crearemos sendos arrays de coordenadas de longitud y latitud con una separación de un grado, expresados en radianes:
lon = np.arange(-180, 181) * np.pi/180
lat = np.arange(-90, 91) * np.pi/180
Lon, Lat = np.meshgrid(lon, lat)
Se necesitará también el radio terrestre expresado en kilómetros, para lo que recurriremos a la librería Astropy:
r = R_earth.to(u.km).value # Radio de la tierra en Km
r
6378.1
Finalmente, las transformaciones trigonométricas:
Xs = r * np.cos(Lon) * np.cos(Lat)
Ys = r * np.sin(Lon) * np.cos(Lat)
Zs = r * np.sin(Lat)
La idea es utilizar la variable z de ETOPO1 para representar las elevaciones mediante un código de colores. Para ello utilizaremos una de las paletas que se han utilizado en el mapa publicado por el NOAA. Las paletas de color utilizadas en ETOPO1 están publicadas aquí. En este caso utilizaremos como colorscale la paleta ETOPO1 CSS3 gradient, que me parece la más fácil de adaptar al formato de las escalas de color de Plotly:
css3g = \
[[0.0, 'rgb( 10, 0,121)'],
[0.0256, 'rgb( 26, 0,137)'],
[0.0513, 'rgb( 38, 0,152)'],
[0.0769, 'rgb( 27, 3,166)'],
[0.1026, 'rgb( 16, 6,180)'],
[0.1282, 'rgb( 5, 9,193)'],
[0.1538, 'rgb( 0, 14,203)'],
[0.1795, 'rgb( 0, 22,210)'],
[0.2051, 'rgb( 0, 30,216)'],
[0.2308, 'rgb( 0, 39,223)'],
[0.2564, 'rgb( 12, 68,231)'],
[0.2821, 'rgb( 26,102,240)'],
[0.3077, 'rgb( 19,117,244)'],
[0.3333, 'rgb( 14,133,249)'],
[0.359, 'rgb( 21,158,252)'],
[0.3846, 'rgb( 30,178,255)'],
[0.4103, 'rgb( 43,186,255)'],
[0.4359, 'rgb( 55,193,255)'],
[0.4615, 'rgb( 65,200,255)'],
[0.4872, 'rgb( 79,210,255)'],
[0.5128, 'rgb( 94,223,255)'],
[0.5385, 'rgb(138,227,255)'],
[0.5641, 'rgb(188,230,255)'],
[0.5641, 'rgb( 51,102, 0)'],
[0.5692, 'rgb( 51,204,102)'],
[0.5744, 'rgb(187,228,146)'],
[0.5897, 'rgb(255,220,185)'],
[0.6154, 'rgb(243,202,137)'],
[0.641, 'rgb(230,184, 88)'],
[0.6667, 'rgb(217,166, 39)'],
[0.6923, 'rgb(168,154, 31)'],
[0.7179, 'rgb(164,144, 25)'],
[0.7436, 'rgb(162,134, 19)'],
[0.7692, 'rgb(159,123, 13)'],
[0.7949, 'rgb(156,113, 7)'],
[0.8205, 'rgb(153,102, 0)'],
[0.8462, 'rgb(162, 89, 89)'],
[0.8718, 'rgb(178,118,118)'],
[0.8974, 'rgb(183,147,147)'],
[0.9231, 'rgb(194,176,176)'],
[0.9487, 'rgb(204,204,204)'],
[0.9744, 'rgb(229,229,229)'],
[1.0, 'rgb(255,255,255)']]
Y finalmente podemos crear la figura Plotly Surface:
data = go.Surface(
x=Xs,
y=Ys,
z=Zs,
colorscale= css3g,
surfacecolor=Z[::-1],
showlegend = False,
showscale=False,
hoverinfo='skip'
)
fig = go.Figure(data)
sin_ejes = dict(showbackground=False,
showgrid=False,
showline=False,
showticklabels=False,
ticks='',
title='',
zeroline=False
)
fig.update_layout(
width=600, height = 600, margin=dict(l=5, r=5, b=5, t=5),
paper_bgcolor = 'black',
scene = dict(
xaxis = sin_ejes,
yaxis = sin_ejes,
zaxis = sin_ejes
)
)
fig.show()
la imagen es interactiva y puede rotarse a voluntad.
Una función muy útil es que puede guardarse la figura interactiva en un fichero:
pio.write_json(fig, 'Datos/Tierra')
Y posteriormente podemos recuperarla y modificarla como vamos a hacer ahora. Esto nos permitirá además reutilizar esta imagen en futuras entradas en este blog.
fig2 = pio.read_json('Datos/Tierra')
A continuación se va a añadir al gráfico anterior la representación gráfica de la órbita de la Estación Espacial Internacional ISS. Seguiremos el método explicado en esta entrada del blog. En primer lugar se creará un objeto EarthSatellite de Skyfield:
ISS_TLE = '''ISS (ZARYA)
1 25544U 98067A 21103.84943184 .00000176 00000-0 11381-4 0 9990
2 25544 51.6434 300.9481 0002858 223.8443 263.8789 15.48881793278621
'''
L0, L1, L2 = ISS_TLE.splitlines()
ISS = EarthSatellite(L1, L2, L0)
Y a continuación creamos el array de posiciones de la ISS en distintos instantes. Pero previamante necesitaremos conocer el periodo del satélite para poder trazar una órbita completa:
rpm = ISS.model.no_kozai # movimiento medio en radianes por minuto
2 * np.pi/rpm #periodo en minutos
92.97029679785253
Es decir, la ISS da una vuelta a La Tierra en aproximadamente 93 minutos.
ts = load.timescale()
tiempos = ts.utc(2021, 4, 13, 0, np.arange(0,94)) # Array de tiempos en minutos
posiciones = ISS.at(tiempos)
pos_ISS = posiciones.position.km
Y para terminar, añadiremos a la figure creada anteriormente una serie de trazas adicionales:
# Añadimos una traza a la figura con la trayectoria de la ISS
fig2.add_trace(go.Scatter3d(x=pos_ISS[0,:], y=pos_ISS[1,:], z=pos_ISS[2,:],
mode='lines', line_color='white', line_width=10,
hoverinfo='skip', showlegend=False))
fig2.update_layout(title_text='Trayectoria de la ISS', title_x=0.5,
title_xanchor='center', title_y=0.9,
title_font_color='white', title_font_size=20)
# Añadimos una traza con el plano ecuatorial
xp =[-1.5 * r, 1.5 * r]
yp = [-1.5 * r, 1.5 * r]
zp = np.zeros((2,2))
color_plano = [[0,"gray"], [1, 'gray']]
fig2.add_trace(go.Surface(x=xp, y=yp, z=zp, colorscale=color_plano,
showscale=False, opacity=0.75, hoverinfo='skip',
showlegend=False))
# Añadimos una traza con el eje terrestre NS
xe = [0,0]
ye = [0,0]
ze = [-1.5 * r, 1.5 * r]
fig2.add_trace(go.Scatter3d(x=xe, y=ye, z=ze,
mode='lines', line_color='red', line_width=5,
hoverinfo='skip', showlegend=False))
fig2.show()
