Justin,
You're in luck - I've had the exact same problem and worked reasonably hard to get an acceptable solution. See the code below... all I've changed is the variable names from to azimuth and elevation to theta and phi... I just had other variables with that naming convention so it all looked nicer.
The trick is grab the highest/lowest elevation point, and copy it across the "top" and "bottom" of the phi/theta space (ie, the poles). You'll see in the first figure which points are copied to where. The "best case scenario" is to have a point ON the poles (in which case no error at all is introduced when that point is copied to various THETA locations), but any point reasonably near to the pole should do just fine.
n=36;
sphSampTh = random('unif',0,2*pi,n,1);
sphSampPhi = random('unif',-pi/2,pi/2,n,1);
r=random('unif',3,7,n,1);
sphSampTh = mod(sphSampTh,2*pi); % Interpret theta from -pi:pi to 0:2pi
% "Wrap" the lateralmost thetas around to the opposite side
wrapSz = pi/4;
wrapThIdxs = find(abs(sphSampTh-pi)>(pi-wrapSz) & abs(sphSampTh-pi)<pi);
wrappedThs = sphSampTh(wrapThIdxs) - 2*pi*sign(sphSampTh(wrapThIdxs)-pi);
% Force the samples taken at the poles to be copied all along the thetaVec but still at the poles.
maxPhiThVec = linspace(0,2*pi,20)';
[~,highPhiIdx] = max(sphSampPhi);
[~,lowPhiIdx] = min(sphSampPhi);
wrapPhiThIdxs = [maxPhiThVec repmat([pi/2 highPhiIdx],size(maxPhiThVec))];
wrapPhiThIdxs = cat(1, wrapPhiThIdxs, [maxPhiThVec repmat([-pi/2 lowPhiIdx],size(maxPhiThVec))]);
% Get a new set of ThPhi combinations that include these additional points, along with indices into
% the sampled coverage map.
gloThPhis = cat(1, [sphSampTh sphSampPhi], [wrappedThs sphSampPhi(wrapThIdxs)], wrapPhiThIdxs(:,1:2));
gloCoverageIdxs = cat(1, (1:n)', wrapThIdxs, wrapPhiThIdxs(:,3));
% NOW we can interp()
DT = DelaunayTri(gloThPhis(:,1),gloThPhis(:,2));
TSinterp = TriScatteredInterp(DT, r(gloCoverageIdxs));
sampThVec = 0:pi/60:2*pi; sampPhiVec = -pi/2:pi/60:pi/2;
[sampThMat, sampPhiMat] = meshgrid(sampThVec, sampPhiVec);
interpdR = TSinterp(sampThMat, sampPhiMat);
% Show the result(s)
figure, plot(sphSampTh,sphSampPhi,'b.',wrappedThs,sphSampPhi(wrapThIdxs),'ro',wrapPhiThIdxs(:,1), wrapPhiThIdxs(:,2),'go')
hold on, rectangle('Position',[0 -.5 2 1]*pi)
title('Scattered spherical sample locations')
xlabel 'Theta', ylabel 'Phi'
figure, imagesc(sampThVec,sampPhiVec,interpdR)
[X Y Z] = sph2cart(sampThMat,sampPhiMat,interpdR);
figure, surf(X,Y,Z)
