function joincst % Generated by BUNDLE: joincst % Created: 1994.10.18.21.03.18 % Copyright (C) 1994 Dr. Charles R. Denham, ZYDECO. % All Rights Reserved. % E-Mail: cdenham@nobska.er.usgs.gov. % WWW: http://priapus.er.usgs.gov. help joincst NL = sprintf('\n'); NEW_FILE = 'fp = fopen([mfunc ''.m''], ''w'');'; ADD_FILE = 'fprintf(fp, mtext);'; END_FILE = 'fclose(fp);'; PROGRESS = 'disp(['' unbundled: '' mfunc])'; mfunc = 'join_cst'; eval(NEW_FILE); mtext = [ ... 'function [new_coast,slen]=join_cst(coast,tol);' NL ... '%% JOIN_CST Makes continuous coastline from fragmented' NL ... '%% segments by joining endpoints of segments less than ' NL ... '%% TOL distance apart. The resulting coastline can be' NL ... '%% used for blanking or filling by adding in a few extra ' NL ... '%% points, and the islands will also be polygons.' NL ... '%% The output coastline is sorted largest' NL ... '%% segment first, followed by sucessively smaller segments.' NL ... '%% ' NL ... '%% USAGE: [new_coast,slen]=join_cst(coast,tol);' NL ... '%%' NL ... '%% Where COAST contains the x and y coastline positions' NL ... '%% in the first two columns, and coastine segments are' NL ... '%% separated by a row containing [NaN NaN] or [-99999. -99999.],' NL ... '%% and TOL is the distance within which coastline segments will' NL ... '%% be joined. SLEN is a vector containing the length of the' NL ... '%% resulting segments. ' NL ... '' NL ... '%%' NL ... '%% Rich Signell | E-mail: rsignell@crusty.er.usgs.gov' NL ... '%% U.S. Geological Survey | Voice : 508-457-2229 ' NL ... '%% 384 Woods Hole Road | FAX : 508-457-2310' NL ... '%% Woods Hole, MA 02543-1598 | WWW : http://crusty.er.usgs.gov/rsignell.html' NL ... '%% ' NL ... '%% requires FIXCOAST.M' NL ... '%%' NL ... 'coast=fixcoast(coast); %%replace -99999. with NaN, eliminate dup NaNs, etc.' NL ... '' NL ... 'plot(coast(:,1),coast(:,2));set(gca,''aspectratio'',[nan 1]);' NL ... 'figure(gcf);' NL ... '' NL ... 'iseg=find(isnan(coast(:,1)));' NL ... ]; eval(ADD_FILE); mtext = [ ... 'nseg=length(iseg)-1;' NL ... 'irem=ones(nseg,1);' NL ... 'z=coast(:,1)+sqrt(-1)*coast(:,2);' NL ... 'zstart=z(iseg(1:nseg)+1);' NL ... 'zstop=z(iseg([2:(nseg+1)])-1);' NL ... 'zends=[zstart(:) zstop(:)];' NL ... '' NL ... 'znew=nan+sqrt(-1)*nan;' NL ... 'k=0;' NL ... '' NL ... '%%while there are remaining segments, process' NL ... 'while(~isempty(find(irem))),' NL ... '%%start at 1st remaining segment' NL ... 'ii=find(irem);' NL ... 'id0=ii(1);' NL ... 'id=id0;' NL ... 'z0=z(iseg(id)+1); %%first point on segment' NL ... 'zind=[(iseg(id)+1):(iseg(id+1)-1)];' NL ... 'line(real(z(zind)),imag(z(zind)),''erasemode'',''none'',''color'',''white'')' NL ... 'zc=z(iseg(id+1)-1); %%last point on segment' NL ... '' NL ... 'irem(id)=0; %%add current segment to list of segments used ' NL ... '' NL ... 'found_next=1;' NL ... 'tried_reverse=0;' NL ... 'while (found_next),' NL ... ' ii=find(irem); %%indices of remaining segments' NL ... ' nrem=length(ii);' NL ... ' dgap=abs((ones(size(zends(ii,:))))*zc-zends(ii,:)); ' NL ... ' id=find(dgapnrem), %%found stop point of segment' NL ... ' id=ii(idt-nrem); %%find which segment we matched ' NL ... ' zi=[(iseg(id+1)-1):-1:(iseg(id)+1)]; %%segment coordinates' NL ... ' else %%found start point of segment' NL ... ' id=ii(idt); %%find which segment we matched ' NL ... ' zi=[(iseg(id)+1):1:(iseg(id+1)-1)]; %%segment coordinates' NL ... ' end' NL ... ' ni=length(zi);' NL ... ' line(real(z(zi)),imag(z(zi)),''erasemode'',''none'',''color'',''white'') %%draw seg' NL ... ' irem(id)=0; %%mask out segment' NL ... ' zc=z(zi(ni)); %%next connecting point' NL ... ' zind=[zind zi]; %%add segment to index list' NL ... ' else' NL ... ' if(~tried_reverse),' NL ... ' zc=z0;' NL ... ' nzind=length(zind);' NL ... ' zind=zind([nzind:-1:1]);' NL ... ' tried_reverse=1;' NL ... ' else' NL ... ' found_next=0;' NL ... ' end' NL ... ' end' NL ... 'end' NL ... 'k=k+1' NL ... 'znew=[znew; nan+sqrt(-1)*nan; z(zind)]; %%add on new contatenated segment ' NL ... 'end' NL ... 'znew=[znew; nan+sqrt(-1)*nan]; %%add nan on end' NL ... 'iseg=find(isnan(znew));' NL ... 'nseg=length(iseg)-1;' NL ... 'slen=(iseg([2:(nseg+1)])-iseg(1:nseg)-2); %%data pts in segments' NL ... '[y,isort]=sort(slen);' NL ... ]; eval(ADD_FILE); mtext = [ ... 'isort=flipud(isort);' NL ... 'zind=[];' NL ... '%%sort largest to smallest' NL ... 'for i=1:nseg;' NL ... ' zind=[zind (iseg(isort(i)):iseg(isort(i)+1))];' NL ... 'end' NL ... 'znew=znew(zind);' NL ... 'new_coast=[real(znew) imag(znew)];' NL ... 'new_coast=fixcoast(new_coast);' NL ... 'iseg=find(isnan(new_coast(:,1)));' NL ... 'nseg=length(iseg)-1;' NL ... 'slen=(iseg([2:(nseg+1)])-iseg(1:nseg)-2); %%data pts in segments' NL ... ]; eval(ADD_FILE); eval(END_FILE); eval(PROGRESS); mfunc = 'fixcoast'; eval(NEW_FILE); mtext = [ ... 'function [coast]=fixcoast(coast)' NL ... '%% function [coast]=fixcoast(coast)' NL ... '%% ' NL ... '%% Fixes coastline is in the format we' NL ... '%% want. Assumes that lon/lat are in the first two columns of' NL ... '%% the matrix coast, and that coastline segments are ' NL ... '%% separated by rows of NaNs (or -99999s). This routine ensures that' NL ... '%% only 1 row of NaNs separates each segment, and makes' NL ... '%% sure that the first and last rows contain NaNs.' NL ... '%%' NL ... 'ind=find(coast==(-99999.));' NL ... 'if(~isempty(ind)),' NL ... ' coast(ind)=coast(ind)*nan;' NL ... 'end' NL ... 'ind=find(isnan(coast(:,1)));' NL ... 'if(~isempty(ind)),' NL ... ' dind=diff(ind);' NL ... ' idup=find(dind==1);' NL ... ' if(~isempty(idup)),' NL ... ' coast(ind(idup),:)=[];' NL ... ' ind(idup)=[];' NL ... ' end' NL ... 'end' NL ... '[m,n]=size(coast);' NL ... 'if(~isnan(coast(1,1))),' NL ... ' coast=[ones(1,n)*nan; coast];' NL ... 'end' NL ... '[m,n]=size(coast);' NL ... 'if(~isnan(coast(m,1))),' NL ... ' coast=[coast; ones(1,n)*nan];' NL ... 'end' NL ... ]; eval(ADD_FILE); eval(END_FILE); eval(PROGRESS);