JOIN_CST Makes continuous coastline from fragmented segments.
USAGE: [new_coast,slen]=join_cst(coast,tol);
Inputs: coast = contains the x and y coastline positions
in the first two columns, and coastine segments are
separated by a row containing [NaN NaN] or [-99999. -99999.],
tol = the distance within which coastline segments will be joined.
Outputs:
new_coast = the new joined coastline, sorted from largest segment
to smallest segment.0001 function [new_coast,slen]=join_cst(coast,tol); 0002 % JOIN_CST Makes continuous coastline from fragmented segments. 0003 % 0004 % USAGE: [new_coast,slen]=join_cst(coast,tol); 0005 % 0006 % Inputs: coast = contains the x and y coastline positions 0007 % in the first two columns, and coastine segments are 0008 % separated by a row containing [NaN NaN] or [-99999. -99999.], 0009 % 0010 % tol = the distance within which coastline segments will be joined. 0011 % 0012 % Outputs: 0013 % new_coast = the new joined coastline, sorted from largest segment 0014 % to smallest segment. 0015 0016 % slen = vector containing the length of the resulting segments. 0017 % 0018 % Description: Makes continous coastline from fragmented segments 0019 % joining endpoints of segments less than 0020 % TOL distance apart. The resulting coastline can be 0021 % used for blanking or filling by adding in a few extra 0022 % points, and the islands will also be polygons. This routine 0023 % draws the new coastline as it works so you can see it's progress. 0024 % 0025 % Need some coastline data? Try out the coastline extractor at: 0026 % 0027 % http://crusty.er.usgs.gov/getcoast.html 0028 % 0029 % Rich Signell | E-mail: rsignell@crusty.er.usgs.gov 0030 % 0031 % requires FIXCOAST.M 0032 % 0033 0034 coast=fixcoast(coast); %replace -99999. with NaN, eliminate dup NaNs, etc. 0035 0036 plot(coast(:,1),coast(:,2));set(gca,'aspectratio',[nan 1]); 0037 figure(gcf); 0038 0039 iseg=find(isnan(coast(:,1))); 0040 nseg=length(iseg)-1; 0041 irem=ones(nseg,1); 0042 z=coast(:,1)+sqrt(-1)*coast(:,2); 0043 zstart=z(iseg(1:nseg)+1); 0044 zstop=z(iseg([2:(nseg+1)])-1); 0045 zends=[zstart(:) zstop(:)]; 0046 0047 znew=nan+sqrt(-1)*nan; 0048 k=0; 0049 0050 %while there are remaining segments, process 0051 while(~isempty(find(irem))), 0052 %start at 1st remaining segment 0053 ii=find(irem); 0054 id0=ii(1); 0055 id=id0; 0056 z0=z(iseg(id)+1); %first point on segment 0057 zind=[(iseg(id)+1):(iseg(id+1)-1)]; 0058 line(real(z(zind)),imag(z(zind)),'erasemode','none','color','white') 0059 zc=z(iseg(id+1)-1); %last point on segment 0060 0061 irem(id)=0; %add current segment to list of segments used 0062 0063 found_next=1; 0064 tried_reverse=0; 0065 while (found_next), 0066 ii=find(irem); %indices of remaining segments 0067 nrem=length(ii); 0068 dgap=abs((ones(size(zends(ii,:))))*zc-zends(ii,:)); 0069 id=find(dgap<tol); 0070 if(~isempty(id)), 0071 idt=find(dgap==min(dgap(:))); % find closest points 0072 idt=idt(1); %take first if equal dist 0073 if(idt>nrem), %found stop point of segment 0074 id=ii(idt-nrem); %find which segment we matched 0075 zi=[(iseg(id+1)-1):-1:(iseg(id)+1)]; %segment coordinates 0076 else %found start point of segment 0077 id=ii(idt); %find which segment we matched 0078 zi=[(iseg(id)+1):1:(iseg(id+1)-1)]; %segment coordinates 0079 end 0080 ni=length(zi); 0081 line(real(z(zi)),imag(z(zi)),'erasemode','none','color','white') %draw seg 0082 irem(id)=0; %mask out segment 0083 zc=z(zi(ni)); %next connecting point 0084 zind=[zind zi]; %add segment to index list 0085 else 0086 if(~tried_reverse), 0087 zc=z0; 0088 nzind=length(zind); 0089 zind=zind([nzind:-1:1]); 0090 tried_reverse=1; 0091 else 0092 found_next=0; 0093 end 0094 end 0095 end 0096 k=k+1 0097 znew=[znew; nan+sqrt(-1)*nan; z(zind)]; %add on new contatenated segment 0098 end 0099 znew=[znew; nan+sqrt(-1)*nan]; %add nan on end 0100 iseg=find(isnan(znew)); 0101 nseg=length(iseg)-1; 0102 slen=(iseg([2:(nseg+1)])-iseg(1:nseg)-2); %data pts in segments 0103 [y,isort]=sort(slen); 0104 isort=flipud(isort); 0105 zind=[]; 0106 %sort largest to smallest 0107 for i=1:nseg; 0108 zind=[zind (iseg(isort(i)):iseg(isort(i)+1))]; 0109 end 0110 znew=znew(zind); 0111 new_coast=[real(znew) imag(znew)]; 0112 new_coast=fixcoast(new_coast); 0113 iseg=find(isnan(new_coast(:,1))); 0114 nseg=length(iseg)-1; 0115 slen=(iseg([2:(nseg+1)])-iseg(1:nseg)-2); %data pts in segments