Index: /issm/trunk-jpl/src/m/classes/basin.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/basin.m	(revision 24001)
+++ /issm/trunk-jpl/src/m/classes/basin.m	(revision 24001)
@@ -0,0 +1,149 @@
+%BASIN class definition
+%
+%   Usage:
+%      basin=basin();
+
+classdef basin
+	properties (SetAccess=public) 
+		boundaries        = {};
+		epsg              = 3426;
+		name              = '';
+	end
+	methods (Static)
+		function self = loadobj(self) % {{{
+			% This function is directly called by matlab when a model object is
+			% loaded. Update old properties here
+		end% }}}
+	end
+	methods
+		function self = basin(varargin) % {{{
+			switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				otherwise
+
+					self=setdefaultparameters(self);
+					options=pairoptions(varargin{:}); 
+			
+					%recover field values: 
+					self.boundaries=getfieldvalue(options,'boundaries',{});
+					self.name=getfieldvalue(options,'name','');
+					self.epsg=getfieldvalue(options,'epsg',3426);
+			end
+		end % }}}
+		function self = setdefaultparameters(self) % {{{
+			self.name='';
+			self.epsg=3426;
+			self.boundaries={};
+
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   basin parameters:'));
+			fielddisplay(self,'name','basin name');
+			fielddisplay(self,'epsg','epsg projection number for the entire basin');
+			fielddisplay(self,'boundaries','list of boundary objects');
+			for i=1:length(self.boundaries),
+				disp(sprintf('             boundary #%i: %s',i,self.boundaries{i}.name));
+			end
+
+		end % }}}
+		function output=outputname(self,varargin) % {{{
+		
+			%recover options
+			options=pairoptions(varargin{:});
+			extension=getfieldvalue(options,'extension',1);
+
+			[path,nme,ext]=fileparts(self.name);
+			if extension,
+				output=[nme ext];
+			else
+				output=nme;
+			end
+		end % }}}
+		function plot(self,varargin) % {{{
+	
+			%add option: 
+			for i=1:length(self.boundaries),
+				self.boundaries{i}.plot('epsg',self.epsg,varargin{:});
+			end
+
+		end % }}}
+		function plot3d(self,varargin) % {{{
+	
+			%add option: 
+			for i=1:length(self.boundaries),
+				self.boundaries{i}.plot3d(varargin{:});
+			end
+
+		end % }}}
+		function out=contour(self,varargin) % {{{
+		
+			%recover options
+			options=pairoptions(varargin{:});
+			x=[]; y=[];
+
+			%go through boundaries, recover edges and project them in the basin epsg reference frame: 
+			for i=1:length(self.boundaries),
+				boundary=self.boundaries{i};
+				contour=boundary.edges();
+				[contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.projection()),sprintf('EPSG:%i',self.epsg));
+				x=[x;contour.x];
+				y=[y;contour.y];
+			end
+			out.x=x;
+			out.y=y;
+		end % }}}
+		function output=shapefilecrop(self,varargin) % {{{
+
+			%recover options
+			options=pairoptions(varargin{:});
+			threshold=getfieldvalue(options,'threshold',.65); %.65 degrees lat,long
+			inshapefile=getfieldvalue(options,'shapefile');
+			outputshapefile=getfieldvalue(options,'outputshapefile','');
+			epsgshapefile=getfieldvalue(options,'epsgshapefile');
+
+			%create list of contours that have critical length > threshold:  (in lat,long)
+			contours=shpread(inshapefile);
+			llist=[];
+			for i=1:length(contours),
+				contour=contours(i);
+				carea=polyarea(contour.x,contour.y);
+				clength=sqrt(carea);
+				if clength<threshold,
+					llist=[llist;i];
+				end
+			end
+			contours(llist)=[];
+
+			%project onto reference frame:
+			if self.epsg~=epsgshapefile,
+				for i=1:length(contours),
+					h=contours(i);
+					[h.x,h.y]=gdaltransform(h.x,h.y,sprintf('EPSG:%i',epsgshapefile),sprintf('EPSG:%i',self.epsg));
+					contours(i).x=h.x;
+					contours(i).y=h.y;
+				end
+			end
+
+			%only keep the contours that are inside the basin (take centroids): 
+			ba=self.contour();
+			flags=zeros(length(contours),1);
+			for i=1:length(contours),
+				h=contours(i); 
+				in=inpolygon(h.x,h.y,ba.x,ba.y); 
+				if ~isempty(find(in==0)),
+					flags(i)=1;
+				end
+			end
+			pos=find(flags);  contours(pos)=[];
+
+			%Two options: 
+			if strcmpi(outputshapefile,''),
+				output=contours;
+			else
+				shpwrite(contours,outputshapefile);
+			end
+
+		end % }}}
+	end
+end
Index: /issm/trunk-jpl/src/m/classes/boundary.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/boundary.m	(revision 24001)
+++ /issm/trunk-jpl/src/m/classes/boundary.m	(revision 24001)
@@ -0,0 +1,216 @@
+%BOUNDARY class definition
+%
+%   Usage:
+%      boundary=boundary();
+
+classdef boundary
+	properties (SetAccess=public) 
+		shppath           = '';
+		shpfilename       = '';
+		orientation       = 'normal';  %other value is 'reverse'
+		epsg              = 4326; %EPSG number, default value is WGS 84 Lat,Long reference frame.
+	end
+	methods (Static)
+		function self = loadobj(self) % {{{
+			% This function is directly called by matlab when a model object is
+			% loaded. Update old properties here
+		end% }}}
+	end
+	methods
+		function self = boundary(varargin) % {{{
+			switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				otherwise
+					self=setdefaultparameters(self);
+					options=pairoptions(varargin{:}); 
+			
+					%recover field values: 
+					self.shppath=getfieldvalue(options,'shppath','');
+					self.shpfilename=getfieldvalue(options,'shpfilename','');
+					self.orientation=getfieldvalue(options,'orientation','normal');
+					self.epsg=getfieldvalue(options,'epsg',4326);
+			end
+		end % }}}
+		function self = setdefaultparameters(self) % {{{
+		self.shppath='';
+		self.shpfilename='';
+		self.orientation='normal';
+		self.epsg=4326;
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   boundary parameters:'));
+
+			fielddisplay(self,'shppath','path to filename for this boundary');
+			fielddisplay(self,'shpfilename','shape file name');
+			fielddisplay(self,'orientation','orientation (default is ''normal'', can be ''reverse'')');
+			fielddisplay(self,'epsg','EPSG number defining projection for the shape file');
+
+		end % }}}
+		function output=name(self) % {{{
+			output=self.shpfilename;
+		end % }}}
+		function output=edges(self) % {{{
+		
+		%read domain:
+		[path,name,ext]=fileparts(self.shpfilename);
+		if ~strcmpi(ext,'shp'),
+			ext='shp';
+		end
+		output=shpread([self.shppath '/' name '.' ext]);
+
+		%do we reverse? 
+		if strcmpi(self.orientation,'reverse'),
+			output.x=flipud(output.x);
+			output.y=flipud(output.y);
+		end
+		end % }}}
+		function output=projection(self) % {{{
+			output=self.epsg;
+		end % }}}
+		function plot(self,varargin) % {{{
+			%recover options
+		
+			options=pairoptions(varargin{:});
+
+			%parse input:
+			figurenumber=getfieldvalue(options,'figure',1);
+			color=getfieldvalue(options,'color','r');
+			linewidth=getfieldvalue(options,'linewidth',1);
+			markersize=getfieldvalue(options,'markersize',1);
+			unitmultiplier=getfieldvalue(options,'unit',1);
+			epsg=getfieldvalue(options,'epsg',4326);
+			radius=getfieldvalue(options,'radius',6371012);
+			aboveground=getfieldvalue(options,'aboveground',1000);
+			offset=getfieldvalue(options,'offset',.1);
+			fontsize=getfieldvalue(options,'fontsize',10);
+
+			%read domain:
+			[path,name,ext]=fileparts(self.shpfilename);
+			if ~strcmpi(ext,'shp'),
+				ext='shp';
+			end
+			domain=shpread([self.shppath '/' name '.' ext]);
+
+			%convert boundary to another reference frame:  {{{
+
+			for i=1:length(domain),
+				try, 
+					[x,y] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),sprintf('EPSG:%i',epsg));
+				catch me
+					disp(me.message());
+					self.disp();
+				end
+				domain(i).x=x; domain(i).y=y;
+			end
+
+			for i=1:length(domain),
+				hold on;
+				if length(x)==1,
+					p=plot(x,y,'k*'); 
+					set(p,'MarkerSize',markersize);
+					t=text(x,y,self.shpfilename,'FontSize',fontsize);
+				else
+					p=plot(x,y,'k-'); 
+					text(sum(x)/length(x),sum(y)/length(y),self.shpfilename,'FontSize',fontsize);
+				end
+				set(p,'Color',color);
+				set(p,'LineWidth',linewidth);
+			end
+			%}}}
+		end % }}}
+		function plot3d(self,varargin) % {{{
+			%recover options
+		
+			options=pairoptions(varargin{:});
+
+			%parse input:
+			figurenumber=getfieldvalue(options,'figure',1);
+			color=getfieldvalue(options,'color','r');
+			linewidth=getfieldvalue(options,'linewidth',1);
+			markersize=getfieldvalue(options,'markersize',1);
+			unitmultiplier=getfieldvalue(options,'unit',1);
+			epsg=getfieldvalue(options,'epsg',4326);
+			radius=getfieldvalue(options,'radius',6371012);
+			aboveground=getfieldvalue(options,'aboveground',1000);
+			offset=getfieldvalue(options,'offset',.1);
+			fontsize=getfieldvalue(options,'fontsize',10);
+
+			%read domain:
+			[path,name,ext]=fileparts(self.shpfilename);
+			if ~strcmpi(ext,'shp'),
+				ext='shp';
+			end
+			domain=shpread([self.shppath '/' name '.' ext]);
+
+			if epsg==4326,
+				%convert boundary to lat,long: {{{
+
+				for i=1:length(domain),
+					try, 
+						[lat,long] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),'EPSG:4326');
+					catch me
+						disp(me.message());
+						self.disp();
+					end
+					domain(i).x=long; domain(i).y=lat;
+				end
+
+				for i=1:length(domain),
+
+					%make sure lat,long are what they are supposed to be: 
+					%if any(domain(i).x>90 | domain(i).x<-90), 
+					%	long=domain(i).x; lat=domain(i).y;
+					%else
+					%	long=domain(i).y; lat=domain(i).x;
+					%end
+
+					%project on x,y,z reference frame.
+					[x,y,z]=AboveGround(domain(i).x,domain(i).y,radius,aboveground);
+					[xt,yt,zt]=AboveGround(domain(i).x+offset,domain(i).y+offset,radius,300000);
+					hold on;
+					if length(x)==1,
+						p=plot3(x,y,z,'k*'); 
+						set(p,'MarkerSize',markersize);
+						t=text(xt,yt,zt,self.shpfilename,'FontSize',fontsize);
+					else
+						p=plot3(x,y,z,'k-'); 
+						mid=floor(length(xt)/2);
+						text(xt(mid),yt(mid),zt(mid),self.shpfilename,'FontSize',fontsize);
+					end
+					set(p,'Color',color);
+					set(p,'LineWidth',linewidth);
+
+				end
+				%}}}
+			else
+				%convert boundary to another reference frame:  {{{
+
+				for i=1:length(domain),
+					try, 
+						[x,y] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),sprintf('EPSG:%i',epsg));
+					catch me
+						disp(me.message());
+						self.disp();
+					end
+					domain(i).x=x; domain(i).y=y;
+				end
+
+				for i=1:length(domain),
+					hold on;
+					if length(x)==1,
+						p=plot(x,y,'k*'); 
+						set(p,'MarkerSize',markersize);
+						t=text(x,y,self.shpfilename,'FontSize',fontsize);
+					else
+						p=plot(x,y,'k-'); 
+						text(sum(x)/length(x),sum(y)/length(y),self.shpfilename,'FontSize',fontsize);
+					end
+					set(p,'Color',color);
+					set(p,'LineWidth',linewidth);
+				end
+				%}}}
+			end
+		end % }}}
+	end
+end