Index: /issm/trunk-jpl/src/m/contrib/larour/exportpoint.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/exportpoint.m (revision 24853) +++ /issm/trunk-jpl/src/m/contrib/larour/exportpoint.m (revision 24853) @@ -0,0 +1,10 @@ +function exportpoint(x,y,shapefilename); + + contours=struct([]); + for i=1:length(x), + contours(i).id=i; + contours(i).Lon=x(i); + contours(i).Lat=y(i); + contours(i).Geometry='Point'; + end + shapewrite(contours,shapefilename); Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 24852) +++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 24853) @@ -16,4 +16,5 @@ shapefileroot = ''; regions = struct(); + o2regions = struct(); element_connectivity = []; glacier_connectivity = []; @@ -32,10 +33,8 @@ 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,5 +45,4 @@ 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,4 +80,11 @@ end %}}} + function units(self) % {{{ + disp('Glacier inventory units: '); + disp(' areas: km^2'); + disp(' mass: Gt'); + + end + %}}} function counter = nglaciers(self,varargin) % {{{ @@ -135,10 +140,10 @@ 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,48 +159,4 @@ 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,108 +229,175 @@ %}}} 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