Im not sure if you need to interpolate at all ... you have a 2d matrix of data with associated vectors for phi and theta? Then you should be able to display it as is using the mapping toolbox.
Lets say you had
theta = [-179:180];
phi = [-89:90];
obs % 360x180 2d array of data
[lat_grid,lon_grid] = meshgrid(theta,phi); % should both be 360x180 to match data
load coast % loads lat and long variables that define the coastline
worldmap('World') % also try axesm as it gives more options
geoshow(lat,long) % draw the coastlines
pcolorm(lat_grid,lon_grid,obs) % a pseudo-color plot of obs data on the projected grid
Hope that helps (alternately, hope you found a solution already)
Cheers, Angus