Strategy question on producing publication plots using matplotlib

Hello,

I'm producing my final publication plots at the moment like that:
- Export a vtk file using

vtk = VTKOutput(ma = self.mesh,coefs=coefs,names=names, filename = directory + filename, subdivision=0, legacy=True) vtk.Do()

- Import that .vtk file again using the vtk library

reader = vtk.vtkUnstructuredGridReader()
        reader.SetFileName(filename)
        reader.ReadAllVectorsOn()
        reader.Update()
    
        data = reader.GetOutput()
    
        # a 2D matrix of size npoints x 3 
        coordinates = np.array(data.GetPoints().GetData())
        # print(coordinates.shape)
        # print(coordinates[0])

        connectivity = np.array(data.GetCells().GetData()).reshape(-1,4)
        # print(connectivity.shape)
        # print(connectivity[0])

        weights = data.GetPointData().GetArray("vec")
        npweights = np.array(weights)

- Do the plotting using matplotlib

fig = plt.figure()
                fig.set_size_inches(10.5, 10.5)
                ax=fig.add_axes([0,0,1,1])
                ax.set_title(label)
    
                tri = mtri.Triangulation(x,y,tri)
    
                ax.triplot(tri, alpha = 0.1)
                ax.scatter(x,y,s=0.5, color='black')
                
                tpc = ax.tripcolor(tri, w, alpha=0.5)
                fig.colorbar(tpc)
                            
                ax.set_aspect('equal')


This method works fine, but seems not super elegant. Are there any suggestions to improve it? Especially, how could I get the  and the  directly from ngsolve fes and Coefficient- / GridFunction?

Thanks!