function latlonoverlay(latlon,numbering)
%LATLONOVERLAY - overlay latitude and longitude lines on current figure
%
%   latstep,lonstep, in latitude and longitude degreees, between two latitudinal, longitudinal profiles.
%   color: [1 1 1] for example
%   resolution: profile resolution ( in lat,lon degrees) 
%   gap: gap (in meters) to plug lat,lon degree numbers;
%
%   Usage:
%      latlonoverlay(latlon,numbering)

%recover latstep and longstep
latstep=latlon{1};
lonstep=latlon{2};
resolution=latlon{3};
color=latlon{4};

if ~iscell(numbering),
	numbering=0;
else
	latgap=numbering{1};
	longap=numbering{2};
	colornumber=numbering{3};
	latangle=numbering{4};
	lonangle=numbering{5};
	numbering=1;
end

%Find lat and lon that fit within the bounds of our image.
lat=-90:resolution:0;
lon=0:resolution:360;

xlimits=xlim;
ylimits=ylim;

for i=1:length(lon),
	if(rem(lon(i),lonstep)==0)
		latitudes=lat';
		longitudes=lon(i)*ones(length(latitudes),1);
		[x,y]=ll2xy(latitudes,longitudes);
		pos=find(x<=xlimits(2) & x>=xlimits(1) & y<=ylimits(2) & y>=ylimits(1));
		if length(pos)<=1,
			continue;
		end
		x=x(pos);y=y(pos);

		l=line(x,y,'Color',color);
		if numbering,
			disp(['Specifiy where to put number for longitude ' num2str(lon(i))]);
			[xcorner,ycorner]=ginput(1);

			%Find nearest point on this profile corresponding to (xcorner,ycorner):
			ind=find_point(x,y,xcorner,ycorner); xcorner=x(ind); ycorner=y(ind);
			if length(x(1:ind))>10,
				xcorner2=x(ind-10); ycorner2=y(ind-10);
			else
				xcorner2=x(ind-1); ycorner2=y(ind-1);
			end
			angle=(180)/pi*atan2((ycorner2-ycorner),(xcorner2-xcorner))+lonangle;
			
			if (xcorner>xlimits(1) & xcorner<xlimits(2) & ycorner>ylimits(1) & ycorner<ylimits(2)),
				th=text(xcorner,ycorner,[num2str(lon(i)) '^{\circ}']);set(th,'Color',colornumber,'Rotation',angle,'FontSize',16,'HorizontalAlignment','center','VerticalAlignment','middle');
				
				%erase line and redraw it in two parts, to leave space for latitude number
				delete(l);
				line(x(1:ind-longap),y(1:ind-longap),'Color',color);hold on;
				line(x(ind+longap:end),y(ind+longap:end),'Color',color);
			end

		end
	end
end

for i=1:length(lat),
	if(rem(lat(i),latstep)==0)
		longitudes=lon';
		latitudes=lat(i)*ones(length(longitudes),1);
		[x,y]=ll2xy(latitudes,longitudes);
		pos=find(x<=xlimits(2) & x>=xlimits(1) & y<=ylimits(2) & y>=ylimits(1));
		if length(pos)<=1,
			continue;
		end
		x=x(pos);y=y(pos);

		l=line(x,y,'Color',color);
		
		if numbering,
			disp(['Specifiy where to put number for latitude ' num2str(lat(i))]);
			[xcorner,ycorner]=ginput(1);

			%Find nearest point on this profile corresponding to (xcorner,ycorner):
			ind=find_point(x,y,xcorner,ycorner); xcorner=x(ind); ycorner=y(ind);
			if length(x(1:ind))>10,
				xcorner2=x(ind-10); ycorner2=y(ind-10);
			else
				xcorner2=x(ind-1); ycorner2=y(ind-1);
			end
			angle=(180)/pi*atan2((ycorner2-ycorner),(xcorner2-xcorner))+latangle;
			
			if (xcorner>xlimits(1) & xcorner<xlimits(2) & ycorner>ylimits(1) & ycorner<ylimits(2)),
				th=text(xcorner,ycorner,[num2str(lat(i)) '^{\circ}']);set(th,'Color',colornumber,'Rotation',angle,'FontSize',16,'HorizontalAlignment','center','VerticalAlignment','middle');
				
				%erase line and redraw it in two parts, to leave space for latitude number
				delete(l);
				line(x(1:ind-latgap),y(1:ind-latgap),'Color',color);hold on;
				line(x(ind+latgap:end),y(ind+latgap:end),'Color',color);
			end

		end

	end
end
