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

../trunk-jpl/src/m/contrib/larour/glacier_inventory.m

@@ -141 +141 @@
                         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 +153 @@
                                 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 % }}}