Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24878)
+++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24879)
@@ -76,4 +76,47 @@
 	end
 	%}}}
+		function readcontours(self) % {{{
+			%we are reading the contours from each shape file for each region, reproject
+			%the latlong in these contours to the local laea projection, and rewrite the shapefile
+			%we a different name extension.
+
+			for i=1:self.nregions,
+				
+				disp(['reading shapefile for region: '  self.regions(i).name]);
+				contours=shaperead([self.root '/' self.regions(i).name '.shp']);
+
+				disp(['concatenating contours (X,Y)']);
+				count=0;
+				for j=1:length(contours),
+					count=count+length(contours(j).X);
+				end
+				xc=zeros(count,1); 
+				yc=zeros(count,1);
+				count=0;
+				for j=1:length(contours),
+					nj=length(contours(j).X);
+					xc((count+j):(count+j+nj-1))=contours(j).X;
+					yc((count+j):(count+j+nj-1))=contours(j).Y;
+					count=count+nj;
+				end
+
+				disp(['projecting (X,Y)']);
+				lat0= fix(mean(self.regions(i).CenLat)); long0=fix(mean(self.regions(i).CenLon));
+				proj=laea(lat0,long0); self.regions(i).proj=proj;
+				[xc,yc]=gdaltransform(xc,yc,'EPSG:4326',proj);
+
+				disp(['plugging back into contours']);
+				count=0;
+				for j=1:length(contours),
+					nj=length(contours(j).X);
+					contours(j).X=xc((count+j):(count+j+nj-1));
+					contours(j).Y=yc((count+j):(count+j+nj-1));
+					count=count+nj;
+				end
+				disp(['saving new contours to disk']);
+				shapewrite(contours,[self.root '/' self.regions(i).name '.laeaproj.shp']);
+			end
+		end
+		%}}}
 		function disp(self) % {{{
 			disp(sprintf('   Glacier inventory:'));