Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24843)
+++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24844)
@@ -16,4 +16,6 @@
 		shapefileroot    = '';
 		regions          = struct();
+		element_connectivity = [];
+		glacier_connectivity = [];
 	end
 	methods
@@ -73,4 +75,48 @@
 		end
 		%}}}
+		function counter = nglaciers(self,varargin) % {{{
+
+			if nargin==1,
+				region=-1; %all regions.
+			else
+				region=varargin{1}; %only one.
+			end
+
+			if region==-1,
+				counter=0;
+				for i=1:self.nregions(),
+					counter=counter+length(self.regions(i).Area);
+				end
+			else
+				counter=length(self.regions(region).Area);
+			end
+		end
+		%}}}
+		function [mpartition,npartition]=partition(self,varargin) % {{{
+			mpartition=zeros(self.nregions(),1);
+			npartition=zeros(self.nregions(),1);
+			counter=0;
+			for i=1:self.nregions(),
+				mpartition(i)=counter+1;
+				npartition(i)=counter+self.nglaciers(i);
+				counter=counter+self.nglaciers(i);
+			end
+
+		end % }}}
+		function [lid,rid]=gidtolid(self,gid) % {{{
+			[mpartition,npartition]=self.partition();
+			for i=1:self.nregions(),
+				if (gid>=mpartition(i) && gid <=npartition(i)),
+					rid=i;
+					lid=gid-mpartition(i)+1;
+					break;
+				end
+			end
+
+		end % }}}
+		function [gid]=lidtogid(self,rid,lid) % {{{
+			[mpartition,npartition]=self.partition();
+			gid=mpartition(rid)+lid-1;
+		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   Glacier inventory:')); 
@@ -79,4 +125,109 @@
 			for i=1:self.nregions(),
 				disp(sprintf('      region %i: ''%s'' %i glaciers ',i,self.regions(i).name,length(self.regions(i).Area)));
+			end
+
+		end % }}}
+		function connectivity_check(self,md) % {{{
+			disp(sprintf('   Checking glacier connectivity:')); 
+			figure(1),clf,hold on;
+			counter=0;
+			
+			%some quick checks:  % {{{
+			if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)),
+				error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!');
+			end
+			if ~isempty(find(self.glacier_connectivity<0)),
+				error('glacier_connectivity table is pointing to elements that are < 0!');
+			end  %}}}
+			
+			for gid=1:length(self.glacier_connectivity),
+				el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:));
+				x1=latel(1); x2=latel(2); x3=latel(3); 
+				y1=longel(1); y2=longel(2); y3=longel(3); 
+
+				[lid,rid]=self.gidtolid(gid);
+				x0=self.regions(rid).CenLat(lid);
+				y0=self.regions(rid).CenLon(lid);
+				
+				if ~isintriangle(x0,x1,x2,x3,y0,y1,y2,y3),
+					%plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-'); text(x1,y1,'1'); text(x2,y2,'2'); text(x3,y3,'3'); plot(x0,y0,'r*');
+					%error(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el));
+					disp(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el));
+					counter=counter+1
+				end
+			end
+
+		end % }}}
+		function connectivity_check2(self,md) % {{{
+			disp(sprintf('   Checking glacier connectivity:')); 
+			figure(1),clf,hold on;
+			
+			%some quick checks:  % {{{
+			if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)),
+				error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!');
+			end
+			if ~isempty(find(self.glacier_connectivity<0)),
+				error('glacier_connectivity table is pointing to elements that are < 0!');
+			end  %}}}
+			
+			for gid=1:length(self.glacier_connectivity),
+				el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:));
+				[lid,rid]=self.gidtolid(gid);
+				lat=self.regions(rid).CenLat(lid);
+				long=self.regions(rid).CenLon(lid);
+				
+				proj=laea(lat,long);
+				[longel;long]
+				[x,y]=gdaltransform([longel;long],[latel;lat],'EPSG:4326',proj);
+				x1=x(1); y1=y(1); x2=x(2); y2=y(2); x3=x(3); y3=y(3); x0=x(4); y0=y(4);
+				
+				%test if x0,y0 is in the 1,2,3 triangle: 
+				epsilon=10^-10;
+				delta=x2*y3-y2*x3-x1*y3+y1*x3+x1*y2-y1*x2;
+
+				%first area coordinate
+				area_1=((y3-y2)*(x0-x3)-(x3-x2)*(y0-y3))./delta
+				%second area coordinate
+				area_2=((x3-x1)*(y0-y3)-(y3-y1)*(x0-x3))./delta
+				%third area coordinate
+				area_3=1-area_1-area_2
+
+				isintriangle=((area_1>=0-epsilon) & (area_2>=0-epsilon) & (area_3>=0-epsilon));
+
+				plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-');
+				plot(x0,y0,'r*');
+				if ~isintriangle,
+					error();
+				end
+
+			end
+
+		end % }}}
+		function connectivity_plot(self,md) % {{{
+			disp(sprintf('   Checking glacier connectivity:')); 
+			figure(1),clf,hold on;
+			
+			%some quick checks:  % {{{
+			if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)),
+				error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!');
+			end
+			if ~isempty(find(self.glacier_connectivity<0)),
+				error('glacier_connectivity table is pointing to elements that are < 0!');
+			end  %}}}
+			
+			for gid=1:length(self.glacier_connectivity),
+				el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:));
+				[lid,rid]=self.gidtolid(gid);
+				lat=self.regions(rid).CenLat(lid); 
+				long=self.regions(rid).CenLon(lid);
+				
+				proj=laea(lat,long);
+				%plot([latel;latel(1)],[longel;longel(1)],'r-*')
+				%plot(lat,long,'k*');
+				[latel;lat]
+				[longel;long]
+				[x,y]=gdaltransform([latel;lat],[longel;long],'EPSG:4326',proj)
+				error;
+
 			end
 
@@ -119,4 +270,112 @@
 		end
 		%}}}
+		function mesh_connectivity(self,mesh) % {{{
+			
+			%initialize glacier_connectivity: 
+			glacier_connectivity=zeros(self.nglaciers,1);
+
+			%assume maximum width for connectivity table and initialize 
+			%as sparse matrix: 
+			ny=3000;
+			element_connectivity=sparse(mesh.numberofelements,ny);
+				
+			disp('Building element and glacier connectivity table: ');
+
+			partition_counter=0;
+			for i=1:self.nregions(),
+				region=self.regions(i);
+
+				disp(['   region ' num2str(i) ': ' region.name]);
+
+				%reproject on a lambert centered on the region: 
+				latel=mesh.lat(mesh.elements)*[1;1;1]/3;
+				longel=mesh.long(mesh.elements)*[1;1;1]/3;
+				for j=1:length(region.CenLat),
+					if mod(j,100)==0,
+						disp(num2str(j/length(region.CenLat)*100));
+					end
+
+					lat0=region.CenLat(j); long0=region.CenLon(j); 
+					%proj=laea(lat0,long0);
+					
+					%to avoid folding of projections from behind the Earth, restrict ourselves
+					%to only locations around a lat,long radius around lat0,long0
+					list_elements=flaglatlongradius(latel,longel,lat0,long0,5);
+
+					%go through these elements and figure out if the glacier is inside: 
+					found=0;
+					for k=1:length(list_elements),
+						el=list_elements(k); 
+						lat=mesh.lat(mesh.elements(el,:));
+						long=mesh.long(mesh.elements(el,:));
+						x0=lat0; x1=lat(1); x2=lat(2); x3=lat(3);
+						y0=long0; y1=long(1); y2=long(2); y3=long(3);
+
+						if isintriangle(x0,x1,x2,x3,y0,y1,y2,y3),
+							%we have found the element where this glacier belongs:
+							found=1;
+							break;
+						end
+					end;
+					if ~found, %try again with a projection and all elements: % {{{
+						proj=laea(lat0,long0);
+						list_elements=(1:mesh.numberofelements)';
+						found=0;
+						[xp,yp]=gdaltransform([mesh.long;long0],[mesh.lat;lat0],'EPSG:4326',proj);
+						x0=xp(end); y0=yp(end);
+						for k=1:length(list_elements),
+							el=list_elements(k); 
+							x=xp(mesh.elements(el,:)); y=yp(mesh.elements(el,:));
+							x1=x(1); x2=x(2); x3=x(3);
+							y1=y(1); y2=y(2); y3=y(3);
+							if isintriangle(x0,x1,x2,x3,y0,y1,y2,y3),
+								%we have found the element where this glacier belongs:
+								found=1;
+								break;
+							end
+						end;
+					end % }}}
+					if ~found, %not found, plot %{{{
+						figure(1),clf,hold on;
+						plot(x0,y0,'r*');
+						for k=1:length(list_elements),
+							el=list_elements(k);
+							x=xp(mesh.elements(el,:)); y=yp(mesh.elements(el,:));
+							x1=x(1); x2=x(2); x3=x(3);
+							y1=y(1); y2=y(2); y3=y(3);
+							plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-');
+						end
+						error('cound not find an element, believe it or not'); %}}}
+					end
+					glacier_connectivity(partition_counter+j)=el;
+				end
+				partition_counter=partition_counter+length(region.CenLat);
+			end
+
+			%Now that we have all the elements for each glacier, we can build a
+			%connectivity table that is the reverse of the glacier connectivity table. 
+			%For each element, this will tell us which glaciers belong to this element. 
+			for j=1:length(glacier_connectivity),
+				el=glacier_connectivity(j);
+				if el==0,
+					continue;
+				end
+				count=element_connectivity(el,ny);
+				if count>ny,
+					error('need to enlarge connectivity table');
+				end
+				element_connectivity(el,count+1)=partition_counter+j;
+				element_connectivity(el,ny)=count+1;
+			end
+
+			%Reduce the number of columns (we did not initially, so we chose an arbitrary ny: 
+			nmax=max(element_connectivity(:,end));
+			element_connectivity=element_connectivity(:,[1:nmax,ny]);
+
+			%Assign pointers: 
+			self.element_connectivity=element_connectivity;
+			self.glacier_connectivity=glacier_connectivity;
+
+		end % }}}
 	end
 end