source: issm/oecreview/Archive/24684-25833/ISSM-24845-24846.diff@ 25834

Index: ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m       (revision 24845)
+++ ../trunk-jpl/src/m/contrib/larour/glacier_inventory.m       (revision 24846)
@@ -141,6 +141,9 @@
                        end  %}}}
                        
                        for gid=1:length(self.glacier_connectivity),
+                               if mod(gid,1000)==0,
+                                       disp(num2str(gid/length(self.glacier_connectivity)*100));
+                               end
                                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); 
@@ -150,56 +153,21 @@
                                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);
+                                       
+                                       %check in local coordinate system: 
+                                       proj=laea(x0,y0);
+                                       [xp,yp]=gdaltransform([y0,y1,y2,y3],[x0,x1,x2,x3],'EPSG:4326',proj);
+                                       x0=xp(1); x1=xp(2); x2=xp(3); x3=xp(4);
+                                       y0=yp(1); y1=yp(2); y2=yp(3); y3=yp(4);
                                
-                               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();
+                                       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
 
                end % }}}