Changeset 24853

Timestamp:

05/13/20 23:17:49 (5 years ago)

Author:

Eric.Larour

Message:

CHG: some serious debugging of the connectivity methods.

Location:

issm/trunk-jpl/src/m/contrib/larour

Files:

• exportpoint.m (added)
• glacier_inventory.m (modified) (7 diffs)

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

@@ -16 +16 @@
                 shapefileroot    = '';
                 regions          = struct();
+                o2regions          = struct();
                 element_connectivity = [];
                 glacier_connectivity = [];
@@ -32 +33 @@
                                 self.shapefileroot=getfieldvalue(options,'shapefileroot');
                                 self.region_names=getfieldvalue(options,'region_names');
-                                regions_zone_shapefile=getfieldvalue(options,'regions_zone_shapefile','');
-
-                                %region region zones first if available: 
-                                if ~strcmpi(regions_zone_shapefile,''),
-                                        region_zones=shpread(regions_zone_shapefile);
-                                end
+                                o2regions_shapefile=getfieldvalue(options,'o2regions_shapefile');
+
+                                %first read o2 regions shapefile: 
+                                self.o2regions=shpread(o2regions_shapefile);
 
                                 %read the shape files and create the regions: 
@@ -46 +45 @@
                                         disp(['reading region: '  self.regions(i).name]);
                                         self.regions(i).id=i;
-                                        self.regions(i).zone=[region_zones(i).x,region_zones(i).y];
                                         contours=shpread([self.shapefileroot '/' self.regions(i).name '.shp']);
                                         
@@ -82 +80 @@
                 end
                 %}}}
+                function units(self) % {{{
+                        disp('Glacier inventory units: ');
+                        disp('   areas: km^2');
+                        disp('   mass: Gt');
+
+                end
+                %}}}
                 function counter = nglaciers(self,varargin) % {{{
 
@@ -135 +140 @@
                 function [lid,rid]=gidtolid(self,gid) % {{{
                         [mpartition,npartition]=self.partition();
+                        lid=zeros(length(gid),1);
+                        rid=zeros(length(gid),1);
                         for i=1:self.nregions(),
-                                if (gid>=mpartition(i) && gid <=npartition(i)),
-                                        rid=i;
-                                        lid=gid-mpartition(i)+1;
-                                        break;
-                                end
+                                pos=find(gid>=mpartition(i) & gid<=npartition(i)); 
+                                rid(pos)=i;
+                                lid(pos)=gid(pos)-mpartition(i)+1;
                         end
 
@@ -154 +159 @@
                         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),
-                                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); 
-
-                                [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),
-                                        
-                                        %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);
-                                
-                                        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
 
@@ -268 +229 @@
                 %}}}
                 function mesh_connectivity(self,mesh) % {{{
-                        
-                        %initialize glacier and element connectivity: 
-                        glacier_connectivity=zeros(self.nglaciers,1);
+
+                        plotr=1;
+
+                        %The way we run this is through the O2 zones defined in o2regions. We go through 
+                        %these  regions, figure out the centroid, figure out how many elements are close to
+                        %this centroid (very important to do this through vertex lat,long, and not through elemnet
+                        % lat,long which can be seriously warped because of the -180 to +180 longitude transition. 
+                        %We then project the region in lamber azimuthal equal area, along with glaciers and elements. 
+                        %Once projected, we figure out which glaciers belong to which element  in the local 
+                        %projection system. 
+
+                        %initialize glacier connectivity: 
+                        self.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);
+                        ny=self.nglaciers(); self.element_connectivity=sparse(mesh.numberofelements,ny);
                                 
                         disp('Building element and glacier connectivity table: ');
 
-                        partition_counter=0;
+                        %O2 regions: 
+                        o1=zeros(self.nglaciers(),1);
+                        o2=zeros(self.nglaciers(),1);
+                        counter=1;
                         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)';
+                                        o1(counter)=region.O1Region(j);
+                                        o2(counter)=region.O2Region(j);
+                                        counter=counter+1;
+                                end
+                        end
+
+                        %Go through O2 regions: 
+                        for i=1:length(self.o2regions),
+                                string=self.o2regions(i).RGI_CODE; 
+                                disp(['progressing with region ' num2str(i) ' ' string]);
+                                offset=findstr(string,'-'); 
+                                o1i=str2num(string(1:offset-1));
+                                o2i=str2num(string(offset+1:end));
+                                glaciers=find(o1==o1i & o2==o2i);
+
+                                if ~isempty(glaciers),
+                                        %find lat,long for laea projection: 
+                                        box=self.o2regions(i).BoundingBox; 
+                                        long0=mean(box(:,1));
+                                        lat0=mean(box(:,2));
+                                        proj=laea(lat0,long0);
+
+                                        %find radius of box: 
+                                        minlat=min(box(:,2)); maxlong=max(box(:,1)); 
+                                        radius=sqrt( (lat0-minlat)^2+(long0-maxlong)^2);
+
+                %some radius adjustment:  {{{
+                switch i,
+                        case 4, radius=40;
+                        case 12, radius=25;
+                        case 33, radius=5;
+                        case 41, radius=75;
+                        case 42, radius=45;
+                        case 68, radius=10;
+                        case 82, radius=30;
+                        otherwise,
+                end % }}}
+
+                                        [lids,rids]=self.gidtolid(glaciers);
+                                        %quick check on the rids, should all be the same number: 
+                                        if min(rids)~=max(rids)error(sprintf('rids should only span on O1 region')); end;
+                                        rid=max(rids);
+
+                                        region=self.regions(rid);
+                                        elements=flaglatlongradiuselements(mesh.elements,mesh.lat,mesh.long,lat0,long0,radius);
+
+                                        if plotr, % {{{
+                                                figure(1),clf; 
+                                                subplot(2,1,1),hold on;
+                                                trisurf(mesh.elements(elements,:),mesh.long,mesh.lat,mesh.lat),view(2); 
+                                                plot3(box(1,1),box(1,2),1000,'r*','MarkerSize',10);
+                                                plot3(box(1,1),box(2,2),1000,'r*','MarkerSize',10);
+
+                                                plot3(box(2,1),box(1,2),1000,'r*','MarkerSize',10);
+                                                plot3(box(2,1),box(2,2),1000,'r*','MarkerSize',10);
+
+                                                plot3(region.CenLon(lids),region.CenLat(lids),1000*ones(length(lids),1),'k*');
+                                        end % }}}
+
+                                        %project lat,long: 
+                                        [x,y]=gdaltransform(mesh.long,mesh.lat,'EPSG:4326',proj);
+                                        [xlid,ylid]=gdaltransform(region.CenLon(lids),region.CenLat(lids),'EPSG:4326',proj);
+
+                                        if plotr, % {{{
+                                                figure(1),subplot(2,1,2), hold on;
+                                                trisurf(mesh.elements(elements,:),x,y,x),view(2); 
+                                                plot3(xlid,ylid,1000*ones(length(lids),1),'k*');
+                                                pause(.1);
+                                        end % }}}
+
+                                        %go through lids: 
+                                        for j=1:length(lids),
                                                 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);
+                                                x0=xlid(j); y0=ylid(j);
+                                                for k=1:length(elements),
+                                                        el=elements(k);
+                                                        x1=x(mesh.elements(el,1)); y1=y(mesh.elements(el,1));
+                                                        x2=x(mesh.elements(el,2)); y2=y(mesh.elements(el,2));
+                                                        x3=x(mesh.elements(el,3)); y3=y(mesh.elements(el,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'); %}}}
+                                                self.glacier_connectivity(glaciers(j))=el;
+                                                if ~found, 
+                                                        error(sprintf('could not find element for glacier %i with lid %i',j,lids(j)));
+                                                end
                                         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);
-                                count=element_connectivity(el,ny);
+                                end
+                        end 
+
+                        %build element connectivity table: 
+                        for j=1:length(self.glacier_connectivity),
+                                el=self.glacier_connectivity(j);
+                                count=self.element_connectivity(el,ny);
                                 if count>ny,
-                                        error('need to enlarge connectivity table');
-                                end
-                                element_connectivity(el,count+1)=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 % }}}
+                                        error('need to enlarge connectivity table to at least');
+                                end
+                                self.element_connectivity(el,count+1)=j;
+                                self.element_connectivity(el,ny)=count+1;
+                        end
+
+                        %Reduce the number of columns (we did not initially, so we chose an arbitrary ny:
+                        nmax=max(self.element_connectivity(:,end));
+                        self.element_connectivity=self.element_connectivity(:,[1:nmax,ny]);
+
+
+                end % }}}
+                function checkconnectivity(self,mesh) % {{{
+
+                        vector=find(self.element_connectivity(:,end));
+
+                        for i=1:length(vector),
+                                el=vector(i);
+
+                                flags=zeros(mesh.numberofelements,1);
+                                flags(el)=1;
+
+                                nglaciers=self.element_connectivity(el,end); 
+                                glaciers=self.element_connectivity(el,1:nglaciers);
+
+                                [lids,rids]=self.gidtolid(glaciers);
+                                lat=zeros(length(glaciers),1);
+                                long=zeros(length(glaciers),1);
+                                for j=1:nglaciers,
+                                        lat(j)=self.regions(rids(j)).CenLat(lids(j));
+                                        long(j)=self.regions(rids(j)).CenLon(lids(j));
+                                end
+
+                                proj=laea(lat(1),long(1));
+
+                                figure(1),clf,hold on;
+                                [x,y]=gdaltransform(mesh.long(mesh.elements(el,:)),mesh.lat(mesh.elements(el,:)),'EPSG:4326',proj);
+                                p=patch('XData',x,'YData',y);  set(p,'FaceColor','red')
+
+                                [x,y]=gdaltransform(long,lat,'EPSG:4326',proj);
+                                plot(x,y,'k*');
+                                pause(.1);
+                        end
+
+                end % }}}
+
         end
 end