Add pfsspy current sheet plotting

[dstansby]

In a pfsspy solution there is a current sheet, which is the surface where $B_{r}$. It would be really nice to be able to plot this in sunkit-pyvista. The tricky bit here will be working out how to extract the surface from the PFSSPy solution.

[wzq215]

Hello, dstansby.
I am also working on visualizing the current sheets, and I think these codes might help.

pfss_in = pfsspy.Input(gong_map, nrho, rss)
pfss_out = pfsspy.pfss(pfss_in)

sg_vect = pfss_out.grid.sg
pg_vect = pfss_out.grid.pg
rg_vect = pfss_out.grid.rg
b_arr = pfss_out.bg

[sg_arr, pg_arr, rg_arr] = np.meshgrid(sg_vect, pg_vect, rg_vect[0:-1])
xg_arr = rg_arr * np.cos(pg_arr) * np.sqrt(1 - sg_arr ** 2)
yg_arr = rg_arr * np.sin(pg_arr) * np.sqrt(1 - sg_arr ** 2)
zg_arr = rg_arr * sg_arr

mesh_g = pyvista.StructuredGrid(xg_arr, yg_arr, zg_arr)
mesh_g.point_data['values'] = b_arr[:, :, 0:-1, 0].ravel(order='F')  # also the active scalars
isos_b0 = mesh_g.contour(isosurfaces=1, rng=[0, 0])
isos_b0.plot(opacity=0.8)

Results:

New to GitHub so I guess it would be better to paste these codes directly rather than risk pulling requests... :)

[dstansby]

😮 that looks amazing, thanks for sharing! I don't have time to help a pull request at the moment, but if you want to have a go at opening a PR we have a guide for newcomers at https://docs.sunpy.org/en/latest/dev_guide/contents/newcomers.html, and maybe someone else would be able to help.

[nabobalis]

If you need some help or advice @wzq215 on how to open a PR, I can provide that.

@dstansby from your point of view, would you like this to be a seperate function we can call to add the current sheet or an addition to the current API?

[wzq215]

Thank you for the guidance! Let me check the guide first.
And there seem to be some problems with the codes above. Hopefully, I can solve that and try opening the PR next week.

[wzq215]

Sorry for getting back so late.
I'm stuck on unifying the coordinates of field lines and the current sheet.
The field lines go across the current sheet in the figure below, which seems not right... I guess something went wrong with the 3d cartesian coordinates for the current sheet (Br).

If you have time please take a look at my code. 👀

# get pfss_out from gong map
gong_map = sunpy.map.Map(gong_fname)
nrho=35
rss=2.5
pfss_in = pfsspy.Input(gong_map, nrho, rss)
pfss_out = pfsspy.pfss(pfss_in)

# extract coordinates of the points & the radial component of B from pfss_out
sg_vect = pfss_out.grid.sg
pg_vect = pfss_out.grid.pg
rg_vect = pfss_out.grid.rg
b_arr = pfss_out.bg
sc_vect = pfss_out.grid.sc
pc_vect = pfss_out.grid.pc
rc_vect = pfss_out.grid.rc
bc_r = pfss_out.bc[0]

# ??? create the 3d cartesian coordinates for Br ???
[s_arr, p_arr, r_arr] = np.meshgrid(sc_vect, pc_vect, rg_vect)
x_arr = np.exp(r_arr) * np.cos(p_arr) * np.sqrt(1 - s_arr ** 2)
y_arr = np.exp(r_arr) * np.sin(p_arr) * np.sqrt(1 - s_arr ** 2)
z_arr = np.exp(r_arr) * s_arr

# plot the isosurface by pyvista
mesh_g = pv.StructuredGrid(x_arr, y_arr, z_arr)
mesh_g.point_data['values'] = bc_r.ravel(order='F')
isos_b0 = mesh_g.contour(isosurfaces=1, rng=[0, 0])

# plot the field lines & the current sheet
tracer = tracing.FortranTracer()
lat = np.deg2rad(np.linspace(-90, 90, 9, endpoint=False))
lon = np.deg2rad(np.linspace(0, 360, 18, endpoint=False))
lat, lon = np.meshgrid(lat, lon, indexing='ij')
lat, lon = lat.ravel() * u.rad, lon.ravel() * u.rad
seeds = SkyCoord(lon, lat, rss * const.R_sun, frame=pfss_out.coordinate_frame)
field_lines = tracer.trace(seeds, pfss_out)

# plot everything
plotter = SunpyPlotter()
plotter.plot_field_lines(field_lines, color_func=my_fline_color_func)
plotter.plotter.add_mesh(pv.Sphere(radius=1))
plotter.plotter.add_mesh(isos_b0, opacity=0.6)
plotter.show()

[nabobalis]

Hi @wzq215, sorry I have not had some time currently to look at your example.

I assume the lines should start at the current sheet and not pass through, that is your issue?

[wzq215]

Hi @nabobalis. Yes, that is the problem.

[TrestanSimon]

Hi @wzq215, apologies for the necropost, but was the issue with this code ever resolved? I am interested in implementing this feature and would be happy to help out.

[wzq215]

Hi @TrestanSimon. Thank you for your attention to this!

In my example, I try to extract the current sheet by extracting the isosurface Br=0 from pfss_output (using mesh.contour() in pyvista).

I think the problem that the field lines pass through the current sheet is caused by:

• the coordinate_frame of pfss_output depends on the input map (either a HeliographicCarrington or HeliographicStonyhurst frame, according to pfsspy Documentation).
• I don't know how to transform the whole pfss_output grid to the coordinate_frame of SunpyPlotter() since pfss_output does not have the property "coord" like the field_line.

Would you happen to know how to transform the coordinate frame of pfss_output?

If I define plotter = SunpyPlotter(coordinate_frame=gong_map.coordinate_frame) and plot the field lines and the current sheet together, the result looks better (figure below). However, some field lines still cross the current sheet (on the left and right limbs). This might suggest that my way of extracting the current sheet is wrong (or incompatible with the fieldline tracer)...

[TrestanSimon]

Hi again @wzq215, I believe I have solved the problem. The pfss_output data just needs to be shifted by the longitude offset of the input map. In FITS files this offset is defined by the keyword CRVAL1.

The code in #53 (comment) would look like

# get pfss_out from gong map
...

# extract coordinates of the points & the radial component of B from pfss_out
...

# get the longitude offset from the input map
lon0 = ((gong_map.meta['crval1'] + 180) * u.deg).to(u.rad).value

# create the 3D Cartesian coordinates for Br
[s_arr, p_arr, r_arr] = np.meshgrid(sc_vect, pc_vect, rg_vect)
x_arr, y_arr, z_arr = pfsspy.coords.strum2cart(r_arr, s_arr, p_arr+lon0)

# plot the isosurface by pyvista
...

# plot the field lines & the current sheet
...

# plot everything
...

where pfsspy.coords.strum2cart was used to convert from the numerical grid coordinates to Cartesian coordinates.

The pfsspy.fieldline class already accounts for this in its coords method by using the property _lon0 of pfss.Output, which is how I figured this out.

[TrestanSimon]

@wzq215 Are you able to open a pull request for this?

[wzq215]

Sorry I've been a little busy this week but I'm working on it :)

[wzq215]

I have opened a pull request (It's my first pull request, and I hope I didn't mess up...)

    # get the offset of the input map (the longitude of the map center)
    lon0 = ((pfss_out.input_map.meta['crval1']) * u.deg).to(u.rad).value
    
    # keep up with self.coordinate_frame
    c = SkyCoord(lon0*u.rad, 0.*u.rad, frame=HeliocentricInertial)
    c = c.transform_to(self.coordinate_frame)
    
    # get grid points. remove the gap
    sc_vect = pfss_out.grid.sc
    pc_vect = np.insert(pfss_out.grid.pc, 0, pfss_out.grid.pc[-1])
    rg_vect = pfss_out.grid.rg
    bc_r = np.insert(pfss_out.bc[0], 0, pfss_out.bc[0][-1], axis=0)
    
    # create 3d grids 
    [s_arr, p_arr, r_arr] = np.meshgrid(sc_vect, pc_vect, rg_vect)
    x_arr, y_arr, z_arr = strum2cart(r_arr, s_arr, p_arr + c.lon.to(u.rad).value)

I'm afraid I'm still not very clear about the coordinates of the pfss output. This code I use works (gets compatible field lines and current sheets in a specified coordinate_frame) but it also seems weird...