Plotting two variables then coloring by a third variable

Question

I have a dataset from an aircraft flight and I am trying to plot the position of the plane (longitude x latitude) then color that line by the altitude of the plan at those coordinates. My code looks like this:

lat_data = np.array( [ 39.916294, 39.87139 , 39.8005  , 39.70801 , 39.64645 , 39.58172 ,
       39.537853, 39.55141 , 39.6787  , 39.796528, 39.91702 , 40.008347,
       40.09513 , 40.144157, 40.090584, 39.96447 , 39.838924, 39.712112,
       39.597103, 39.488377, 39.499096, 39.99354 , 40.112175, 39.77281 ,
       39.641186, 39.51512 , 39.538853, 39.882736, 39.90413 , 39.811333,
       39.73279 , 39.65676 , 39.584026, 39.5484  , 39.54484 , 39.629486,
       39.96    , 40.07143 , 40.187405, 40.304718, 40.423153, 40.549305,
       40.673313, 40.794548, 40.74402 , 40.755558, 40.770306, 40.73574 ,
       40.795086, 40.774628] )

long_data = np.array( [ -105.13034 , -105.144104, -105.01132 , -104.92708 , -104.78505 ,
       -104.6449  , -104.49255 , -104.36578 , -104.32623 , -104.31285 ,
       -104.32199 , -104.41774 , -104.527435, -104.673935, -104.81152 ,
       -104.82184 , -104.81882 , -104.81314 , -104.74657 , -104.78108 ,
       -104.93442 , -104.98039 , -105.0168  , -105.04967 , -105.056564,
       -105.03639 , -105.13429 , -105.05214 , -105.17435 , -105.070526,
       -104.93587 , -104.80029 , -104.65973 , -104.50339 , -104.33972 ,
       -104.21634 , -103.96216 , -103.84808 , -103.72534 , -103.60455 ,
       -103.48926 , -103.376495, -103.25937 , -103.10858 , -103.08469 ,
       -103.24878 , -103.4169  , -103.53073 , -103.23694 , -103.41254 ] )

altitude_data = np.array( [1.6957603e+00,  1.9788861e+00,  1.8547169e+00,  1.8768315e+00,
        1.9633590e+00,  2.0504241e+00,  2.1115899e+00,  2.1085002e+00,
        1.8621666e+00,  1.8893014e+00,  1.8268168e+00,  1.7574688e+00,
        1.7666028e+00,  1.7682364e+00,  1.8120643e+00,  1.7637002e+00,
        1.8054264e+00,  1.9149075e+00,  2.0173934e+00,  2.0875392e+00,
        2.1486480e+00,  1.8622510e+00,  1.7937366e+00,  1.8748144e+00,
        1.9063262e+00,  1.9397615e+00,  2.1261981e+00,  2.0180094e+00,
        1.9827688e+00, -9.9999990e+06,  1.8933343e+00,  1.9615903e+00,
        2.1000245e+00,  2.1989927e+00,  2.3200927e+00, -9.9999990e+06,
        4.0542388e+00,  4.0591464e+00,  4.0597038e+00,  4.3395977e+00,
        4.6702847e+00,  5.0433373e+00,  5.2824092e+00,  5.2813010e+00,
        5.2735353e+00,  5.2784677e+00,  5.2784038e+00,  5.2795196e+00,
        4.9482727e+00,  4.2531524e+00] )

import matplotlib as plt    

fig, ax1 = plt.subplots( figsize = ( 10, 10 ) )
ax1.plot( long_data, lat_data, alpha = .4)
ax1.scatter( long_data, lat_data, c = altitude_data )
plt.show()

Which gives us this track: .

Is there a way to consolidate the data into one line that plots the location of the aircraft and adjusts the color for the elevation?

While plotting a line and a scatter together works, it does not look very good when I put in all the data (n = 2400 ). Thanks!

Answer 1

It looks like if you want to use a Line2D object, you're stuck with a single color per object. As a workaround, you could plot each line segment as a set of (first order linearly) interpolated segments and color each of those by its corresponding infinitesimal value.

It looks like this functionality is contained in a LineCollection instance, however I just went for a more quick and dirty approach below.

For extra credit, since we're talking about geospatial data here, why not use cartopy to plot your data? That way you can have a "basemap" which gives you some reference. After all, if it's worth plotting, it's worth plotting beautifully.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import cartopy
import cartopy.crs as ccrs

import numpy as np
import scipy
from scipy import interpolate

import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt

### clean data
filter_inds   = np.where(np.abs(altitude_data) < 100)
lat_data      = lat_data[filter_inds]
long_data     = long_data[filter_inds]
altitude_data = altitude_data[filter_inds]

# =============== plot

plt.close('all')
plt.style.use('dark_background') ## 'default'
fig = plt.figure(figsize=(1500/100, 1000/100))
#ax1 = plt.gca()

lon_center = np.mean(long_data); lat_center = np.mean(lat_data)

ax1 = plt.axes(projection=ccrs.Orthographic(central_longitude=lon_center, central_latitude=lat_center))
ax1.set_aspect('equal')

scale = 3 ### 'zoom' with smaller numbers
ax1.set_extent((lon_center-((0.9*scale)), lon_center+((0.7*scale)), lat_center-(0.5*scale), lat_center+(0.5*scale)), crs=ccrs.PlateCarree())

### states
ax1.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural', scale='10m', facecolor='none', name='admin_1_states_provinces_shp'), zorder=2, linewidth=1.0, edgecolor='w')

ax1.add_feature(cartopy.feature.RIVERS.with_scale('10m'), zorder=2, linewidth=1.0, edgecolor='lightblue')
ax1.add_feature(cartopy.feature.LAKES.with_scale('10m'), zorder=2, linewidth=1.0, edgecolor='gray')

### download counties from https://prd-tnm.s3.amazonaws.com/StagedProducts/Small-scale/data/Boundaries/countyl010g_shp_nt00964.tar.gz
### untar with : tar -xzf countyl010g_shp_nt00964.tar.gz

try:
    reader = cartopy.io.shapereader.Reader('countyl010g.shp')
    counties = list(reader.geometries())
    COUNTIES = cartopy.feature.ShapelyFeature(counties, ccrs.PlateCarree())
    ax1.add_feature(COUNTIES, facecolor='none', alpha=0.5, zorder=2, edgecolor='gray')
except:
    pass

#norm = matplotlib.colors.Normalize(vmin=altitude_data.min(), vmax=altitude_data.max())
norm = matplotlib.colors.Normalize(vmin=1.0, vmax=6.0)
cmap = matplotlib.cm.viridis
mappableCmap = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)

# ===== plot line segments individually for gradient effect

for i in range(long_data.size-1):
    long_data_this_segment = long_data[i:i+2]
    lat_data_this_segment  = lat_data[i:i+2]
    altitude_data_this_segment  = altitude_data[i:i+2]
    
    ### create linear interp objects
    ### scipy doesnt like when the data isn't ascending (hence the flip)
    
    try:
        spl_lon = scipy.interpolate.splrep(altitude_data_this_segment, long_data_this_segment, k=1)
        spl_lat = scipy.interpolate.splrep(altitude_data_this_segment, lat_data_this_segment,  k=1)
    except:
        long_data_this_segment = np.flip(long_data_this_segment)
        lat_data_this_segment = np.flip(lat_data_this_segment)
        altitude_data_this_segment = np.flip(altitude_data_this_segment)
        spl_lon = scipy.interpolate.splrep(altitude_data_this_segment, long_data_this_segment, k=1)
        spl_lat = scipy.interpolate.splrep(altitude_data_this_segment, lat_data_this_segment,  k=1)
    
    ### linearly resample on each segment
    nrsmpl=100
    altitude_data_this_segment_rsmpl = np.linspace(altitude_data_this_segment[0],altitude_data_this_segment[1],nrsmpl)
    long_data_this_segment_rsmpl = scipy.interpolate.splev(altitude_data_this_segment_rsmpl, spl_lon)
    lat_data_this_segment_rsmpl = scipy.interpolate.splev(altitude_data_this_segment_rsmpl, spl_lat)
    
    for j in range(long_data_this_segment_rsmpl.size-1):
        
        long_data_this_segment_2 = long_data_this_segment_rsmpl[j:j+2]
        lat_data_this_segment_2  = lat_data_this_segment_rsmpl[j:j+2]
        altitude_data_this_segment_2  = altitude_data_this_segment_rsmpl[j:j+2]
        
        ax1.plot(long_data_this_segment_2, lat_data_this_segment_2, transform=ccrs.PlateCarree(), c=mappableCmap.to_rgba(np.mean(altitude_data_this_segment_2)), zorder=3, linestyle='solid', alpha=0.8, lw=5.0)

# =====

### plot the actual data points as a scatter plot
pts = ax1.scatter(long_data, lat_data, transform=ccrs.PlateCarree(), alpha=1.0, marker='o', c=mappableCmap.to_rgba(altitude_data), edgecolor='w', zorder=4)

cbar = fig.colorbar(mappable=mappableCmap, ax=ax1, orientation='vertical', fraction=0.046, pad=0.04)
cbar.set_label(r'$Altitude$ [units]', fontsize=20)
cbar.ax.tick_params(labelsize=16)
cbar.set_ticks(np.linspace(1.0, 6.0, 5+1), update_ticks=True)
cbar.set_ticklabels([ ('%0.1f' % x) for x in cbar.get_ticks() ])

fig.tight_layout()
fig.savefig('flightPath.png',dpi=100)
plt.show()