diff --git a/enlilviz/enlil.py b/enlilviz/enlil.py index 6114551..ecdbd61 100644 --- a/enlilviz/enlil.py +++ b/enlilviz/enlil.py @@ -106,6 +106,31 @@ def get_satellite_data(self, satellite, var, coord=None): varname += '_' + coord return self.ds.loc[{'satellite': satellite}][varname] + @_validate_satellite + def get_satellite_fieldline(self, satellite, time): + """Returns the position of the fieldline passing through the satellite. + + Parameters + ---------- + satellite : str + Satellite of interest. + time : datetime-like + Time of interest. Looks for nearest time. + + Returns + ------- + r : float + Radial positions of the fieldline (AU) + lat : float + Latitudinal positions of the fieldline (deg). + lon : float + Longitudinal positions of the fieldline (deg). + """ + # Fieldlines are index based on 't' not 'earth_t' + ds = self.ds.sel({'t': time}, method='nearest') + ds = ds.loc[{'satellite': satellite}] + return (ds['fld_pos_r'], ds['fld_pos_lat'], ds['fld_pos_lon']) + def get_slice(self, var, slice_plane, time=None): """Get a 2D slice of data. diff --git a/enlilviz/io.py b/enlilviz/io.py index 642b174..ffb0a3b 100644 --- a/enlilviz/io.py +++ b/enlilviz/io.py @@ -250,13 +250,13 @@ def _transform_variable(da): # Position validation is coordinate dependent so do those # conversions first - if da.name == 'pos_r': + if 'pos_r' in da.name: da /= _unit_conversion['AU'] da.attrs['units'] = 'AU' - elif da.name == 'pos_lat': + elif 'pos_lat' in da.name: da.values = np.rad2deg(np.pi/2 - da) da.attrs['units'] = 'degrees_north' - elif da.name == 'pos_lon': + elif 'pos_lon' in da.name: da.values = np.rad2deg(da - np.pi) da.attrs['units'] = 'degrees_east' diff --git a/enlilviz/plotting/plots.py b/enlilviz/plotting/plots.py index d24685f..ec00eeb 100644 --- a/enlilviz/plotting/plots.py +++ b/enlilviz/plotting/plots.py @@ -222,6 +222,18 @@ def _init_plot(self): color=SAT_COLORS[sat], markeredgecolor='k', markersize=10, zorder=2) self.plot_data[sat] = marker + # Fieldlines + pos = run.get_satellite_fieldline(sat, self.time) + mask = np.logical_or(pos[0].values <= RMIN, + pos[0].values >= RMAX) + + lon_field = np.ma.masked_array(np.deg2rad(pos[2]), mask=mask) + r_field = np.ma.masked_array(pos[0], mask=mask) + + line, = ax.plot(lon_field, r_field, + color=SAT_COLORS[sat], + linestyle='-.', zorder=1) + self.plot_data[sat + '_fieldline'] = line # Circle on Earth (1AU) # ax.axvline(0., c='k', zorder=1) @@ -258,6 +270,17 @@ def update(self): pos = run.get_satellite_position(sat, self.time) self.plot_data[sat].set_data(np.deg2rad(pos[2]), pos[0]) + # Fieldlines + pos = run.get_satellite_fieldline(sat, self.time) + mask = np.logical_or(pos[0].values <= RMIN, + pos[0].values >= RMAX) + + lon_field = np.ma.masked_array(np.deg2rad(pos[2]), mask=mask) + r_field = np.ma.masked_array(pos[0], mask=mask) + + self.plot_data[sat + '_fieldline'].set_data( + lon_field, r_field) + class LongitudeSlice(_BasePlot): """A longitude polar slice plot, which is in latitude-radial coordinates. @@ -310,6 +333,19 @@ def _init_plot(self): markeredgecolor='k', markersize=10, zorder=2) self.plot_data[sat] = marker + # Fieldlines + pos = run.get_satellite_fieldline(sat, self.time) + mask = np.logical_or(pos[0].values <= RMIN, + pos[0].values >= RMAX) + + lat_field = np.ma.masked_array(np.deg2rad(pos[1]), mask=mask) + r_field = np.ma.masked_array(pos[0], mask=mask) + + line, = ax.plot(lat_field, r_field, + color=SAT_COLORS[sat], + linestyle='-.', zorder=1) + self.plot_data[sat + '_fieldline'] = line + # x == theta, y == r circle_points = np.linspace(0, 2*np.pi) ax.plot(circle_points, np.ones(len(circle_points)), c='k', @@ -342,6 +378,16 @@ def update(self): for sat in ['Earth']: pos = run.get_satellite_position(sat, self.time) self.plot_data[sat].set_data(np.deg2rad(pos[1]), pos[0]) + # Fieldlines + pos = run.get_satellite_fieldline(sat, self.time) + mask = np.logical_or(pos[0].values <= RMIN, + pos[0].values >= RMAX) + + lat_field = np.ma.masked_array(np.deg2rad(pos[1]), mask=mask) + r_field = np.ma.masked_array(pos[0], mask=mask) + + self.plot_data[sat + '_fieldline'].set_data( + lat_field, r_field) class RadialSlice(_BasePlot):