Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
286 lines (253 sloc) 10.5 KB
function varargout=ellesmereg(res,buf)
% XY=ELLESMEREG(res,buf)
%
% Finds the coordinates of Ellesmere Island, taking into account
% the fact that it is next to Greenland and Baffin Island, which
% causes problems when you buffer.
%
% INPUT:
%
% res 0 The standard, default values
% N Splined values at N times the resolution
% buf Distance in degrees that the region outline will be enlarged
% by BUFFERM, not necessarily integer, possibly negative
% [default: 0]
%
% OUTPUT:
%
% XY Closed-curved coordinates that result from this process
%
% NOTES: We start with a region outline created by encircling glaciers in
% the Randolph Glacier Inventory 3.2 (http://www.glims.org/RGI/). This
% outline is buffered as requested. Since buffered regions infringe on
% neighboring regions, we subtract the outlines for Greenland and Baffin
% Island. This can cause some artifacts near some of the boundary
% intersections, so this function employs the same brushing technique as the
% basin functions for Antarctica. We make the region, plot it, remove
% erroneous points using the brushing tool, and save this new outline.
%
% Last modified by charig at princeton.edu, 09/04/2014
% Last modified by fjsimons-at-alum.mit.edu, 09/23/2014
defval('res',0)
defval('buf',0)
if ~isstr(res) % Not a demo
% The directory where you keep the base coordinates
whereitsat=fullfile(getenv('IFILES'),'GLACIERS','RGI_3_2');
% Revert to original name if unbuffered
if res==0 && buf==0
fnpl=fullfile(whereitsat,'Ellesmere.mat');
elseif buf==0;
fnpl=fullfile(whereitsat,sprintf('%s-%i.mat','Ellesmere',res));
elseif buf~=0
fnpl=fullfile(whereitsat,sprintf('%s-%i-%g.mat','Ellesmere',res,buf));
end
% If you already have a file
if exist(fnpl,'file')==2
load(fnpl)
if nargout==0
plot(XY(:,1),XY(:,2),'k-'); axis equal; grid on
else
varns={XY};
varargout=varns(1:nargout);
end
else
% You are about to make a file
% Do we buffer? Not here, so do it regular
if buf==0
if res==0
% First part
load(fullfile(getenv('IFILES'),'Glaciers','RGI_3_2','Ellesmere.mat'));
% We've already checked this file when we made it that it is correct
else
XY=ellesmere(0);
XYb=bezier(XY,res);
XY=XYb;
end
end
if buf ~=0
% Check if we have this buffer already but only at base resolution, and then change
% the res on that file
fnpl2=fullfile(whereitsat,sprintf('%s-%i-%g.mat','Ellesmere',0,buf));
if exist(fnpl2,'file')==2
load(fnpl2)
XYb=bezier(XY,res);
XY=XYb;
else
% We make a new buffer
disp('Buffering the coastlines... this may take a while');
XY=ellesmere(10);
if buf > 0
inout='out';
else
inout='in';
end
%%%
% First we make the thing next door (Greenland)
%%%
XY1 = ellesmere(10);
if buf>0.7
XY1 = ellesmere(10,buf-0.5);
LonB=XY1(:,1); LatB=XY1(:,2);
else
% Here we use buf=0.2 as a minumum
[LatB,LonB] = bufferm(XY1(:,2),XY1(:,1),0.2,'out');
% Note that, if due to BEZIER there might be a pinched-off loop in
% the XY you will get an extra NaN and will need to watch it
% If this shouldn't work, keep it all unchanged in other words
try
% You'll need a line for every possible version behavior
% Note how POLY2CW has disappeared from BUFFERM
if sign(buf)<0 || ~isempty(strfind(version,'2010a'))
% Take the last bit of non-NaNs; there might have been pinched
% off sections in the original
LonB=LonB(indeks(find(isnan(LonB)),'end')+1:end);
LatB=LatB(indeks(find(isnan(LatB)),'end')+1:end);
elseif ~isempty(strfind(version,'2011a')) || ~isempty(strfind(version,'2012a'))
LonB=LonB(1:find(isnan(LonB))-1);
LatB=LatB(1:find(isnan(LatB))-1);
end
catch
disp('BUFFERM failed to buffer as expected')
end
end % end buf > 0.7
% Now subtract a smaller Ellesmere from Greenland
% Note: Say you ask for buf=1.0. Then we just made Greenland1.0
% and we subtract from it Ellesmere0.5. This is considered the
% thing next door. If we just subtracted a Greenland1.0, we
% would lose large sections of Ellesmere.
LonB(LonB<0) = LonB(LonB<0)+360;
XY1 = [LonB LatB];
XY2 = greenland(10,buf);
[x1,y1] = polybool('subtraction',XY2(:,1),XY2(:,2),XY1(:,1),XY1(:,2));
% Do you need to see it to understand? Use the code below
%path(path,'~/src/m_map');
%figure
%m_proj('oblique mercator','longitudes',[318 318],'latitudes',[90 50],'aspect',1.0);
%m_grid;
%m_coast('color','k');
%m_line(x1(:,1),y1(:,2),'color','magenta','linestyle','-');
%%%
% Now we get Ellesmere
%%%
[LatB,LonB] = bufferm(XY(:,2),XY(:,1),abs(buf),inout);
% Check for weird loops due to bezier
try
% You'll need a line for every possible version behavior
% Note how POLY2CW has disappeared from BUFFERM
if sign(buf)<0 || ~isempty(strfind(version,'2010a'))
% Take the last bit of non-NaNs; there might have been pinched
% off sections in the original
LonB=LonB(indeks(find(isnan(LonB)),'end')+1:end);
LatB=LatB(indeks(find(isnan(LatB)),'end')+1:end);
elseif ~isempty(strfind(version,'2011a')) || ~isempty(strfind(version,'2012a'))
LonB=LonB(1:find(isnan(LonB))-1);
LatB=LatB(1:find(isnan(LatB))-1);
end
catch
disp('BUFFERM failed to buffer as expected')
end
% Periodize our way
LonB(LonB<0) = LonB(LonB<0)+360;
% Now subtract the thing next door we made above so we don't
% overlap Greenland
[x,y] = polybool('subtraction',LonB,LatB,x1,y1);
%%%
% The part for Baffin
%%%
% Also subtract the Baffin island part
% The function for Baffin has internal code that prevents it from
% overlapping Ellesmere, so we need not worry about that here.
if buf>0.7
XY3 = baffing(10,buf-0.5);
[x,y] = polybool('subtraction',x,y,XY3(:,1),XY3(:,2));
else
XY3 = baffing(10,0.2);
[x,y] = polybool('subtraction',x,y,XY3(:,1),XY3(:,2));
end
%%%
% Assembly complete!
%%%
% Now we look at the new piece, and we know we must fix some edges
% where the things intersected.
hdl1=figure;
plot(x,y);
title('This plot is used to edit the coastlines.')
fprintf(['The functions ELLESMERE has paused, and made a plot \n'...
' of the current coastlines. These should have some artifacts that \n'...
'you want to remove. Here are the instructions to do that:\n \n'])
fprintf(['DIRECTIONS: Select the data points you want to remove with \n'...
'the brush tool. Then right click and remove them. After you have\n'...
' finished removing the points you want, select the entire curve \n'...
'with the brush tool, and type return. The program will save the \n'...
'currently brushed data in a variable, and then make another plot \n'...
'for you to confirm you did it right.\n'])
keyboard
% Get the brushed data from the plot
pause(0.1);
hBrushLine = findall(hdl1,'tag','Brushing');
brushedData = get(hBrushLine, {'Xdata','Ydata'});
brushedIdx = ~isnan(brushedData{1});
brushedXData = brushedData{1}(brushedIdx);
brushedYData = brushedData{2}(brushedIdx);
figure
plot(brushedXData,brushedYData)
title('This figure confirms the new data you selected with the brush.')
fprintf(['The newest figure shows the data you selected with the brush \n'...
'tool after you finished editing. If this is correct, type return.\n'...
' If this is incorrect, type dbquit and run this program again to redo.\n'])
keyboard
XY = [brushedXData' brushedYData'];
end % end if fnpl2 exist
end % end if buf>0
% Save the file
save(fnpl,'XY')
% Output
varns={XY};
varargout=varns(1:nargout);
end % end if fnpl exist
elseif strcmp(res,'demomake')
% This demo illustrates the proper order that is needed to make these
% coordinate files, since there are dependancies.
% The coordinates need to be created in order of increasing buffer,
% and ellesmere before baffin
XY = ellesmere(10,0.2);
XY = baffin(10,0.2);
XY = ellesmere(10,0.5);
XY = baffin(10,0.5);
XY = ellesmere(10,1);
XY = baffin(10,1);
XY = ellesmere(10,1.5);
XY = baffin(10,1.5);
XY = ellesmere(10,2.0);
XY = baffin(10,2.0);
elseif strcmp(res,'demo1')
path(path,'~/src/m_map');
XY1 = ellesmere(10);
XY2 = ellesmere(10,0.2);
XY3 = greenland(1,0.5);
figure
m_proj('oblique mercator','longitudes',[318 318],'latitudes',[90 50],'aspect',1.0);
m_grid;
m_coast('color','k');
% Original
m_line(XY1(:,1),XY1(:,2),'color','magenta','linestyle','-');
% Buffered
m_line(XY2(:,1),XY2(:,2),'color','blue','linestyle','-');
m_line(XY3(:,1),XY3(:,2),'color','green','linestyle','-');
elseif strcmp(res,'demo2')
path(path,'~/src/m_map');
XY1 = ellesmere(10,0.2);
XY2 = ellesmere(10,0.5);
XY3 = greenland(1,0.5);
[x1,y1] = polybool('subtraction',XY3(:,1),XY3(:,2),XY1(:,1),XY1(:,2));
[x2,y2] = polybool('subtraction',XY2(:,1),XY2(:,2),x1,y1);
figure
m_proj('oblique mercator','longitudes',[318 318],'latitudes',[90 50],'aspect',1.0);
m_grid;
m_coast('color','k');
% Original
m_line(x2,y2,'color','magenta','linestyle','-');
% Buffered
%m_line(XY2(:,1),XY2(:,2),'color','blue','linestyle','-');
%m_line(XY3(:,1),XY3(:,2),'color','green','linestyle','-');
end
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.