mu_coast.m

A mapping package for Matlab:这是一款功能十分强大的地理绘图工具包

  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;

edg=9999+zeros(length(x),1);
tol=1e-10;

switch MAP_VAR_LIST.rectbox,
  case {'on','off'},
    i=abs(x-xl(1))<tol;
    edg(i)=(y(i)-yl(1))/diff(yl);

    i=abs(x-xl(2))<tol;
    edg(i)=2+(yl(2)-y(i))/diff(yl);

    i=abs(y-yl(1))<tol;
    edg(i)=3+(xl(2)-x(i))/diff(xl);

    i=abs(y-yl(2))<tol;
    edg(i)=1+(x(i)-xl(1))/diff(xl);

  case 'circle',
    i=abs(x.^2 + y.^2 - rl2)<tol;
    edg(i)=-atan2(y(i),x(i));   % use -1*angle so that numeric values
                                % increase in CW direction
end;


%%
function [ncst,k,Area]=get_coasts(optn,file);
%
%  GET_COASTS  Loads various GSHHS coastline databases and does some preliminary
%              processing to get things into the form desired by the patch-filling
%              algorithm.
%
% Changes; 3/Sep/98 - RP: better decimation, fixed bug in limit checking.

global MAP_PROJECTION MAP_VAR_LIST

llim=rem(MAP_VAR_LIST.longs(1)+360,360)*1e6;
rlim=rem(MAP_VAR_LIST.longs(2)+360,360)*1e6;
tlim=MAP_VAR_LIST.lats(2)*1e6;
blim=MAP_VAR_LIST.lats(1)*1e6;

mrlim=rem(MAP_VAR_LIST.longs(2)+360+180,360)-180;
mllim=rem(MAP_VAR_LIST.longs(1)+360+180,360)-180;
mtlim=MAP_VAR_LIST.lats(2);
mblim=MAP_VAR_LIST.lats(1);

% decfac is for decimation of areas outside the lat/long bdys.
% Sizes updated for gshhs v1.3
switch optn(1),
  case 'f',   % 'full' (undecimated) database
    ncst=NaN+zeros(2063513,2);Area=zeros(153545,1);k=ones(153546,1);
    decfac=12500;
  case 'h',
    ncst=NaN+zeros(2063513,2);Area=zeros(153545,1);k=ones(153546,1);
    decfac=2500;
  case 'i',
    ncst=NaN+zeros(493096,2);Area=zeros(41529,1);k=ones(41530,1);
    decfac=500;
  case 'l',
    ncst=NaN+zeros(101207,2);Area=zeros(10775,1);k=ones(10776,1);
    decfac=100;
  case 'c',
    ncst=NaN+zeros(14884,2);Area=zeros(1870,1);k=ones(1871,1);
    decfac=20;
end;
fid=fopen(file,'r','ieee-be');

if fid==-1,
  warning(sprintf(['Coastline file ' file ...
          ' not found \n(Have you installed it? See the M_Map User''s Guide for details)' ...
          '\n ---Using default coastline instead']));
  load m_coasts
  return
end;


Area2=Area;

% Read the File header
%%[A,cnt]=fread(fid,9,'int32');
 [A,cnt]=get_gheader(fid);



l=0;
while cnt>0,

 % A: 1:ID, 2:num points, 3:land/lake etc., 4-7:w/e/s/n, 8:area, 9:greenwich crossed.
 
 C=fread(fid,A(2)*2,'int32'); % Read all points in the current segment.

 %For Antarctica the lime limits are 0 to 360 (exactly), thus c==0 and the
 %line is not chosen for (e.g. a conic projection of part of Antarctica)
 % Fix 30may/02
 if A(5)==360e6, A(5)=A(5)-1; end;
 
 a=rlim>llim;  % Map limits cross longitude jump? (a==1 is no)
 b=A(9)<65536; % Cross boundary? (b==1 if no).
 c=llim<rem(A(5)+360e6,360e6); 
 d=rlim>rem(A(4)+360e6,360e6);
 e=tlim>A(6) & blim<A(7);
 
 % This test checks whether the lat/long box containing the line overlaps that of
 % the map. There are various cases to consider, depending on whether map limits
 % and/or the line limits cross the longitude jump or not.
 
%% if e & (  a&( b&c&d | ~b&(c|d)) | ~a&(~b | (b&(c|d))) ),
 if e & (  (a&( (b&c&d) | (~b&(c|d)) )) | (~a&(~b | (b&(c|d))) ) ),
 
   l=l+1;
 
   x=C(1:2:end)*1e-6;y=C(2:2:end)*1e-6;
 %%  plot(x,y);pause;
   
   %  make things continuous (join edges that cut across 0-meridian)

   dx=diff(x);
%%fprintf('%f %f\n', max(dx),min(dx))
%   if A(9)>65536 | any(dx)>356 | any(dx<356),
%  if ~b | any(dx>356) | any(dx<-356),
     x=x-360*cumsum([x(1)>180;(dx>356) - (dx<-356)]);
%   end;

   % Antarctic is a special case - extend contour to make nice closed polygon
   % that doesn't surround the pole.   
   if abs(x(1))<1 & abs(y(1)+68.9)<1,
     y=[-89.9;-78.4;y(x<=-180);y(x>-180);   -78.4;-89.9*ones(18,1)];
     x=[  180; 180 ;x(x<=-180)+360;x(x>-180);-180; [-180:20:160]'];
   end;

   % First and last point should be the same.
   
   if x(end)~=x(1) | y(end)~=y(1), x=[x;x(1)];y=[y;y(1)]; end;

   % get correct curve orientation for patch-fill algorithm.
   
   Area2(l)=sum( diff(x).*(y(1:(end-1))+y(2:end))/2 );
   Area(l)=A(8)/10;

   if rem(A(3),2)==0; 
     Area(l)=-abs(Area(l)); 
     if Area2(l)>0, x=x(end:-1:1);y=y(end:-1:1); end;
   else
     if Area2(l)<0, x=x(end:-1:1);y=y(end:-1:1); end; 
   end;

   % Here we try to reduce the number of points.
   
   xflag=0;
   if max(x)>180, % First, save original curve for later if we anticipate
     sx=x;sy=y;   % a 180-problem.
     xflag=1;
   end;
   
   % Look for points outside the lat/long boundaries, and then decimate them
   % by a factor of about 'decfac' (don't get rid of them completely because that
   % can sometimes cause problems when polygon edges cross curved map edges).
   
   tol=.2;   
  
   % Do y limits, then x so we can keep corner points.
   
   nn=(y>mtlim+tol) | (y<mblim-tol);
     % keep one extra point when crossing limits, also the beginning/end point.
   nn=logical(nn-min(1,([0;diff(nn)]>0)+([diff(nn);0]<0)));
   nn([1 end])=0;
     % decimate vigorously
   nn=nn & rem(1:length(nn),decfac)'~=0;
   x(nn)=[];y(nn)=[];
         
   if mrlim>mllim,  % no wraparound
       % sections of line outside lat/long limits
     nn=(x>mrlim+tol | x<mllim-tol) & y<mtlim & y>mblim;
    else            % wraparound case
     nn=(x>mrlim+tol & x<mllim-tol ) & y<mtlim & y>mblim;
   end;
   nn=logical(nn-min(1,([0;diff(nn)]>0)+([diff(nn);0]<0)));nn([1 end])=0;
   nn=nn & rem(1:length(nn),decfac)'~=0;
   x(nn)=[];y(nn)=[];
   
   % Move all points "near" to map boundaries.
   % I'm not sure about the wisdom of this - it might be better to clip
   % to the boundaries instead of moving. Hmmm. 
      
   y(y>mtlim+tol)=mtlim+tol;
   y(y<mblim-tol)=mblim-tol;
   if mrlim>mllim,   % Only clip long bdys if I can tell I'm on the right
                     % or left (i.e. not in wraparound case)
     x(x>mrlim+tol)=mrlim+tol;
     x(x<mllim-tol)=mllim-tol;
   end;   
   
 %% plot(x,y);pause;
 
   k(l+1)=k(l)+length(x)+1;
   ncst(k(l)+1:k(l+1)-1,:)=[x,y];
   ncst(k(l+1),:)=[NaN NaN];
   
   % This is a little tricky...the filling algorithm expects data to be in the
   % range -180 to 180 deg long. However, there are some land parts that cut across
   % this divide so they appear at +190 but not -170. This causes problems later...
   % so as a kludge I replicate some of the problematic features at 190-360=-170.
   % Small islands are just duplicated, for the Eurasian landmass I just clip
   % off the eastern part.
   
   if xflag,
     l=l+1;Area(l)=Area(l-1); 
     if abs(Area(l))>1e5,
       nn=find(sx>180);nn=[nn;nn(1)];
       k(l+1)=k(l)+length(nn)+1;
       ncst(k(l)+1:k(l+1)-1,:)=[sx(nn)-360,sy(nn)];
     else   % repeat the island at the other edge.
       k(l+1)=k(l)+length(sx)+1;
       ncst(k(l)+1:k(l+1)-1,:)=[sx-360,sy];
     end;
     ncst(k(l+1),:)=[NaN NaN];
   end;
 end;
 
 
 %%[A,cnt]=fread(fid,9,'int32');
 [A,cnt]=get_gheader(fid);
 
end;

fclose(fid);

%%plot(ncst(:,1),ncst(:,2));pause;clf;
%size(ncst)
%size(Area)
%size(k)


ncst((k(l+1)+1):end,:)=[];  % get rid of unused part of data matrices
Area((l+1):end)=[];
k((l+2):end)=[];
 
%%%
function [A,cnt]=get_gheader(fid);
% Reads the gshhs file header
% 
% A bit of code added because header format changed with version 1.3.

% This works for version 1.2, but not 1.3.

[A,cnt]=fread(fid,9,'int32');

if cnt==9 & A(9)==3,  % we have version 1.3, this would be one of 0,1,65535,65536 in v 1.2
  % Read one more byte
  A2=fread(fid,1,'int32');
  % This is the easiest way not to break existing code.
  A(9)=A2;
end;