mu_coast.m

m_map是加拿大学者编写的一个在matlab上绘地图的软件

function [ncst,Area,k]=mu_coast(optn,varargin);
% MU_COAST Add a coastline to a given map.
%         MU_COAST draw a coastline as either filled patches (slow) or
%         lines (fast) on a given projection. It uses a coastline database with
%         a resolution of about 1/4 degree. 
%
%         It is not meant to be called directly by the user, instead it contains
%         code common to various functions in the m_map directory.
%
%         Calling sequence is MU_COAST(option,arguments) where
%
%         Option string  
%           c,l,i,h,f :  Accesses various GSHHS databases. Next argument is
%                        GSHHS filename.
%
%           u[ser] :   Accesses user-specified coastline file (a mat-file of
%                      data saved from a previous MU_COAST call). Next argument
%                      is filename
%
%           v[ector] : Uses input vector of data. Next argument is the 
%                     data in the form of a nx2 matrix of [longitude latitude].
%                     Patches must have first and last points the same. In a vector,
%                     different patches can be separated by NaN.
%
%           d[efault] :  Accesses default coastline.
%
%         The arguments given above are then followed by (optional) arguments 
%         specifying lines or patches, in the form:
%
%          optional arguments:  <line property/value pairs>, or
%                               'line',<line property/value pairs>.
%                               'patch',<patch property/value pairs>.
%                               'speckle',width,density,<line property/value pairs>.
%
%         If no or one output arguments are specified then the coastline is drawn, with
%         patch handles returned. 
%         If 3 output arguments are specified in the calling sequence, then no drawing
%         is done. This can be used to save subsampled coastlines for future use
%         with the 'u' option, for fast drawing of multiple instances of a particular
%         coastal region.
%
%
%    
%         See also M_PROJ, M_GRID     

% Rich Pawlowicz ([e-mail]) 2/Apr/1997
%
% Notes: 15/June/98 - fixed some obscure problems that occured sometimes
%                     when using conic projections with large extents that
%                     crossed the 180-deg longitude line.
%        31/Aug/98  - added "f" gshhs support (Thanks to RAMO)
%        17/June/99 - 'line' option as per manual (thanks to Brian Farrelly)
%        3/Aug/01   - kludge fix for lines to South Pole in Antarctica (response
%                     to Matt King).
%        30/May/02  - fix to get gshhs to work in Antarctica.
%        15/Dec/05  - speckle additions
%        21/Mar/06  - handling of gshhs v1.3 (developed from suggestions by
%                     Martin Borgh)
%        26/Nov/07 - changed 'finite' to 'isfinite; after warnings
%

% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.


global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION),
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end;


% m_coasts.mat contains 3 variables:
% ncst: a Nx2 matrix of [LONG LAT] line segments, each of which form
%       a closed contour, separated by NaN
%  k=[find(isnan(ncst(:,1)))];
% Area: a vector giving the areas of the different patches. Both ncst
%     and Area should be ordered with biggest regions first. Area can
%     be computed as follows:
%
%  Area=zeros(length(k)-1,1);
%  for i=1:length(k)-1,
%    x=ncst([k(i)+1:(k(i+1)-1) k(i)+1],1);
%    y=ncst([k(i)+1:(k(i+1)-1) k(i)+1],2);
%    nl=length(x);
%    Area(i)=sum( diff(x).*(y(1:nl-1)+y(2:nl))/2 );
%  end;
%
%     Area should be >0 for land, and <0 for lakes and inland seas.

switch optn(1),
  case {'c','l','i','h','f'},  
    [ncst,k,Area]=get_coasts(optn,varargin{1});
    varargin(1)=[];
  case  'u',                   
    eval(['load ' varargin{1} ' -mat']);
    varargin(1)=[];
  case 'v',
    ncst=[NaN NaN;varargin{1};NaN NaN];
    varargin(1)=[];
    k=[find(isnan(ncst(:,1)))];  % Get k
    Area=ones(size(k));          % Make dummy Area vector (all positive).
  otherwise
    load m_coasts
end;

% If all we wanted was to extract a sub-coastline, return.
if nargout==3,
  return;
end;

% Handle wrap-arounds (not needed for azimuthal and oblique projections)

switch MAP_PROJECTION.routine,
 case {'mp_cyl','mp_conic','mp_tmerc'}
  if MAP_VAR_LIST.longs(2)<-180,
   ncst(:,1)=ncst(:,1)-360;
  elseif MAP_VAR_LIST.longs(1)>180,
   ncst(:,1)=ncst(:,1)+360;
  elseif MAP_VAR_LIST.longs(1)<-180,
   Area=[Area;Area];
   k=[k;k(2:end)+k(end)-1];
   ncst=[ncst;[ncst(2:end,1)-360 ncst(2:end,2)]];
   % This is kinda kludgey - but sometimes adding all these extra points
   % causes wrap-around in the conic projection, so we want to limit the
   % longitudes to the range needed. However, we don't just clip them to
   % min long because that can cause problems in trying to decide which way
   % curves are oriented when doing the fill algorithm below. So instead
   % I sort of crunch the scale, preserving topology.
   %
   % 12/Sep/2006 - in the gsshs_crude database we have a lot of long skinny
   % things which interact badly with this - so I offset the scrunch 2 degrees
   % away from the bdy
   nn=ncst(:,1)<MAP_VAR_LIST.longs(1)-2;
   ncst(nn,1)=(ncst(nn,1)-MAP_VAR_LIST.longs(1)+2)/100+MAP_VAR_LIST.longs(1)-2;
  elseif MAP_VAR_LIST.longs(2)>180,
   Area=[Area;Area];
   k=[k;k(2:end)+k(end)-1];
   ncst=[ncst;[ncst(2:end,1)+360 ncst(2:end,2)]];
   % Ditto.
   nn=ncst(:,1)>MAP_VAR_LIST.longs(2)+2;
   ncst(nn,1)=(ncst(nn,1)-MAP_VAR_LIST.longs(2)-2)/100+MAP_VAR_LIST.longs(2)+2;
  end;
end;


if length(varargin)>0,
  if strcmp(varargin(1),'patch'),
    optn='patch';
  end
  if strcmp(varargin(1),'speckle'),
    optn='speckle';
  end
  if strcmp(varargin(1),'line'),
    optn='line';
    varargin=varargin(2:end); % ensure 'line' does not get passed to line
  end
else
  optn='line';
end;



switch optn,
 case {'patch','speckle'}

  switch MAP_VAR_LIST.rectbox,
    case 'on',
      xl=MAP_VAR_LIST.xlims;
      yl=MAP_VAR_LIST.ylims;
      [X,Y]=m_ll2xy(ncst(:,1),ncst(:,2),'clip','on');
      oncearound=4;
    case 'off',
      xl=MAP_VAR_LIST.longs;
      yl=MAP_VAR_LIST.lats;
      X=ncst(:,1);
      Y=ncst(:,2);
      [X,Y]=mu_util('clip','on',X,xl(1),X<xl(1),Y);
      [X,Y]=mu_util('clip','on',X,xl(2),X>xl(2),Y);
      [Y,X]=mu_util('clip','on',Y,yl(1),Y<yl(1),X);
      [Y,X]=mu_util('clip','on',Y,yl(2),Y>yl(2),X);
      oncearound=4;
    case 'circle',
      rl=MAP_VAR_LIST.rhomax;
      [X,Y]=m_ll2xy(ncst(:,1),ncst(:,2),'clip','on');
      oncearound=2*pi;    
  end;

  p_hand=zeros(length(k)-1,1); % Patch handles
  
  for i=1:length(k)-1,
    x=X(k(i)+1:k(i+1)-1);
    fk=isfinite(x);
    if any(fk),
      y=Y(k(i)+1:k(i+1)-1);
%% if i>921, disp('pause 1'); pause; end; 
      nx=length(x);
      if Area(i)<0, x=flipud(x);y=flipud(y); fk=flipud(fk); end;
%clf
%line(x,y,'color','m');
      st=find(diff(fk)==1)+1;
      ed=find(diff(fk)==-1);
%length(x),
%ed,
%st
       if length(st)<length(ed) | isempty(st), st=[ 1;st]; end;
       if length(ed)<length(st),               ed=[ed;nx]; end;
%ed
%st
      if  ed(1)<st(1),
        if c_edge(x(1),y(1))==9999,
          x=x([(ed(1)+1:end) (1:ed(1))]);
          y=y([(ed(1)+1:end) (1:ed(1))]);
          fk=isfinite(x);
          st=find(diff(fk)==1)+1;
          ed=[find(diff(fk)==-1);nx];
          if length(st)<length(ed), st=[1;st]; end
        else
          ed=[ed;nx];
          st=[1;st];
        end;
      end;
%ed
%st
      % Get rid of 2-point curves (since often these are out-of-range lines with
      % both endpoints outside the region of interest)
      k2=(ed-st)<3;
      ed(k2)=[]; st(k2)=[];
%%[XX,YY]=m_ll2xy(x(st),y(st),'clip','off');
%line(x,y,'color','r','linest','-');
%line(x(st),y(st),'marker','o','color','r','linest','none');
%line(x(ed),y(ed),'marker','o','color','g','linest','none');
      edge1=c_edge(x(st),y(st));
      edge2=c_edge(x(ed),y(ed));
%-edge1*180/pi
%-edge2*180/pi
      mi=1;
      while length(st)>0,
        if mi==1,
          xx=x(st(1):ed(1));
          yy=y(st(1):ed(1));
        end;
        estart=edge2(1);
        s_edge=edge1;
        s_edge(s_edge<estart)=s_edge(s_edge<estart)+oncearound;
%s_edge,estart
        [md,mi]=min(s_edge-estart);
        switch MAP_VAR_LIST.rectbox,
          case {'on','off'},
            for e=floor(estart):floor(s_edge(mi)),
              if e==floor(s_edge(mi)), xe=x(st([mi mi])); ye=y(st([mi mi])); 
              else  xe=xl; ye=yl; end;
              switch rem(e,4),
                case 0,
                  xx=[xx; xx(end*ones(10,1))];
                  yy=[yy; yy(end)+(ye(2)-yy(end))*[1:10]'/10 ];
                case 1,
                  xx=[xx; xx(end)+(xe(2)-xx(end))*[1:10]'/10 ];
                  yy=[yy; yy(end*ones(10,1))];
                case 2,
                  xx=[xx; xx(end*ones(10,1))];
                  yy=[yy; yy(end)+(ye(1)-yy(end))*[1:10]'/10;];
                case 3,
                  xx=[xx; xx(end)+(xe(1)-xx(end))*[1:10]'/10 ];
                  yy=[yy; yy(end*ones(10,1))];
              end;
            end;
          case 'circle',
            if estart~=9999,
%s_edge(mi),estart
              xx=[xx; rl*cos(-[(estart:.1:s_edge(mi))]')];
              yy=[yy; rl*sin(-[(estart:.1:s_edge(mi))]')];
            end;
        end;
        if mi==1, % joined back to start
	  if strcmp(optn,'patch'),
            if strcmp(MAP_VAR_LIST.rectbox,'off'), 
              [xx,yy]=m_ll2xy(xx,yy,'clip','off'); 
            end;
             if Area(i)<0,
              p_hand(i)=patch(xx,yy,varargin{2:end},'facecolor',get(gca,'color'));
            else
              p_hand(i)=patch(xx,yy,varargin{2:end});
            end;
	  else  % speckle
            if ~strcmp(MAP_VAR_LIST.rectbox,'off'),  % If we were clipping in
	       [xx,yy]=m_xy2ll(xx,yy);                % screen coords go back
	    end;  
            if Area(i)>0,
               p_hand(i)=m_hatch(xx,yy,'speckle',varargin{2:end});
	    else   
	       p_hand(i)=m_hatch(xx,yy,'outspeckle',varargin{2:end});
            end;
	  end;  
%%%if i>921, disp(['paused-2 ' int2str(i)]);pause; end;
          ed(1)=[];st(1)=[];edge2(1)=[];edge1(1)=[];
        else
          xx=[xx;x(st(mi):ed(mi))];
          yy=[yy;y(st(mi):ed(mi))];
          ed(1)=ed(mi);st(mi)=[];ed(mi)=[];
          edge2(1)=edge2(mi);edge2(mi)=[];edge1(mi)=[];
        end;
%%if i>921, disp(['paused-2 ' int2str(i)]);pause; end;
      end;
    end; 
  end;  

 otherwise,
 
  % This handles the odd points required at the south pole by any Antarctic
  % coastline by setting them to NaN (for lines only)
  ii=ncst(:,2)<=-89.9;
  if any(ii), ncst(ii,:)=NaN; end;
  
  [X,Y]=m_ll2xy(ncst(:,1),ncst(:,2),'clip','on');
 
  % Get rid of 2-point lines (these are probably clipped lines spanning the window)
  fk=isfinite(X);        
  st=find(diff(fk)==1)+1;
  ed=find(diff(fk)==-1);
  k=find((ed-st)==1);
  X(st(k))=NaN;

  p_hand=line(X,Y,varargin{:}); 

end;

ncst=p_hand;

%-----------------------------------------------------------------------
function edg=c_edge(x,y);
% C_EDGE tests if a point is on the edge or not. If it is, it is
%        assigned a value representing it's position oon the perimeter
%        in the clockwise direction. For x/y or lat/long boxes, these
%        values are
%           0 -> 1 on left edge
%           1 -> 2 on top
%           2 -> 3 on right edge
%           3 -> 4 on bottom
%        For circular boxes, these values are the -ve of the angle
%        from center.

global MAP_VAR_LIST

switch MAP_VAR_LIST.rectbox,
  case 'on',
    xl=MAP_VAR_LIST.xlims;
    yl=MAP_VAR_LIST.ylims;
  case 'off',
    xl=MAP_VAR_LIST.longs;
    yl=MAP_VAR_LIST.lats;
  case 'circle',
    rl2=MAP_VAR_LIST.rhomax^2;
end;