Layout algorithms in tools like networkx, igraph, and gephi will associate coordinates with your nodes which you should be able to access fairly easily. Once you have those coordinates, you just need to plot your pie-charts on top of the relevant node location. Alternative, these tools also support using external images as node markers, so instead of building the plots in the same script you could build the pie charts separately, save them to disk, and then associate them with nodes when you draw the graph.
I've never seen an "out-of-the-box" solution for this specific kind of graphic, but it shouldn't be too hard to do this yourself. You just need to figure out how to access the layout coordinates. If you clarify what your preferred analytic environment is and/or graph analysis tool, I can give you more specific advice.
EDIT: I managed to find the code that was used to build the chart in that paper. I searched the paper for "we used" and found this in the acknowledgements:
We are especially indebted to Aaron Clauset and James Fowler for thorough readings of a draft of this manuscript and to Christina Frost for developing some of the graph visualizations we used.
Searching "Christina Frost UNC" led me to this page which contains a collection of graph visualization tools for matlab. The one you are looking for is at the bottom:
drawForceCPie.m. The site is super slow, but it eventually shared the code with me. Here it is for posterity in case the site crashes:
map=[.7 .7 .7; map];
% for i=1:length(colorsu);
idcolors = map(ceil(Rcolor),:);
% for i=1:length(colorsu),
% if colorsu(i)==1
set(h,'color',[.5 .5 .5]*(1-str(i,j)));
nodes_percom = length(find(gnu(i)==gn));
points = 40;
x = pos(1);
y = pos(2);
last_t = 0;
for i = 1:length(percents)
end_t = last_t + percents(i)*points;
tlist = [last_t ceil(last_t):floor(end_t) end_t];
xlist = [0 (radius*cos(tlist*2*pi/points)) 0] + x;
ylist = [0 (radius*sin(tlist*2*pi/points)) 0] + y;
last_t = end_t;
tlist = [0:points];
xlist = x+radius*cos(tlist*2*pi/points);
ylist = y+radius*sin(tlist*2*pi/points);
function mat = commAdjMatrixSparse(groups, A)
% Creates a community adjacency matrix using
% groups from the output for reccurrcommsNew2Sparse, A is the adjacency matrix
% 0's on the diagonal, other elements consist of the total number of
% connections between the two communities
% Last Modified by ALT 20 June 2007
rows = max(communities);
for i = 1:rows
for j = 1:rows
if(i ~= j)
comm1 = find(communities==i);
comm2 = find(communities==j);
mat(j, i) = sum(sum(A(comm1, comm2)));
function [communities cut] = findcommunitiesatcut(groups,cut)
% Gives the community numbers at a requested cut or level in the groups vector,
% if the cut number is not valid the program changes it to a valid one.
% Uses a groups vector and a scalar cut number, gives communities and the cut number,
% which is needed when cut is changed.
%Last modified by ALT, 20 June 2007
disp(['That is too many cuts! I have changed the cut number.']);
disp(['Negative numbers dont work, I have changed the cut number to the max'])
%Identify distinct group values and number of cut levels in dendrogram:
cutdiff=diffnumbers(cuts+1-cut); %NOTE THERE IS NO ERROR CHECKING HERE, ASSUMED VALID CUT NUMBER
%Define communities by replacing the groupnumbers values in groups with the
%corresponding commnumbers values, component by component.
%Is there a more efficient way to specify this in MATLAB?
If you use this code for research, I believe this is the citation you should reference (in addition to the UNC webpage that hosted the code):
"Visualization of communities in networks," Amanda L. Traud, Christina Frost, Peter J. Mucha, and Mason A. Porter, Chaos 19, 041104 (2009).