Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 24853) +++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (revision 24854) @@ -6,15 +6,20 @@ % where varargin is a variable list of options: % +% Usage: +% rgi = glacier_inventory('root',shapefileroot,... +% 'filenames',region_names,... +% 'boxes', region_boxes); % Example: -% rgi = glacier_inventory('shapefileroot',shapefileroot,'region_names',region_names); +% rgi = glacier_inventory('root','~/ModelData',... +% 'filenames',{'01_rgi60_Alaska','02_rgi60_WesternCanadaUS'},... +% 'boxes','00_rgi60_O2Regions.shp'); % -% +% Watch out boxes do not have to match the region shapefiles. Example: O1regions vs O2 regions shapefiles. classdef glacier_inventory < handle properties (SetAccess=public) %Model fields - region_names = {}; - shapefileroot = ''; + root = ''; regions = struct(); - o2regions = struct(); + boxes = struct(); element_connectivity = []; glacier_connectivity = []; @@ -24,133 +29,50 @@ function self = glacier_inventory(varargin) % {{{ - if nargin==0, - self=setdefaultparameters(self); - else - self=setdefaultparameters(self); - - options=pairoptions(varargin{:}); - - self.shapefileroot=getfieldvalue(options,'shapefileroot'); - self.region_names=getfieldvalue(options,'region_names'); - o2regions_shapefile=getfieldvalue(options,'o2regions_shapefile'); - - %first read o2 regions shapefile: - self.o2regions=shpread(o2regions_shapefile); - - %read the shape files and create the regions: - counter=0; - self.regions=struct(); - for i=1:self.nregions(), - self.regions(i).name=self.region_names{i}; - disp(['reading region: ' self.regions(i).name]); - self.regions(i).id=i; - contours=shpread([self.shapefileroot '/' self.regions(i).name '.shp']); - - %we can't keep all the info: build arrays of centroids: - O1Region=zeros(length(contours),1); - O2Region=zeros(length(contours),1); - Area=zeros(length(contours),1); - CenLon=zeros(length(contours),1); - CenLat=zeros(length(contours),1); - for j=1:length(contours), - O1Region(j)=str2num(contours(j).O1Region); - O2Region(j)=str2num(contours(j).O2Region); - Area(j)=contours(j).Area; - CenLon(j)=contours(j).CenLon; - CenLat(j)=contours(j).CenLat; - end - - self.regions(i).Area=Area; - self.regions(i).O1Region=O1Region; - self.regions(i).O2Region=O2Region; - self.regions(i).CenLat=CenLat; - self.regions(i).CenLon=CenLon; - self.regions(i).lids=[1:length(contours)]'; - self.regions(i).gids=self.regions(i).lids+counter; - counter=counter+length(contours); - end - end - end - %}}} - function inv = setdefaultparameters(inv) % {{{ - end - %}}} - function n = nregions(self) % {{{ - n=length(self.region_names); - end - %}}} - function units(self) % {{{ - disp('Glacier inventory units: '); - disp(' areas: km^2'); - disp(' mass: Gt'); - - end - %}}} - function counter = nglaciers(self,varargin) % {{{ - - if nargin==1, - region=-1; %all regions. - else - region=varargin{1}; %only one. - end - - if region==-1, - counter=0; - for i=1:self.nregions(), - counter=counter+length(self.regions(i).Area); - end - else - counter=length(self.regions(region).Area); - end - end - %}}} - function name=ridtoname(self,rid) % {{{ - - fullname=self.region_names{rid}; - name=fullname(10:end); - - end % }}} - function totalarea=area(self,varargin) % {{{ - region=-1; - totalarea=0; - if nargin==2, - region=varargin{1}; - end - if region==-1, - %figure out the areas of everybody: - for i=1:self.nregions(), - totalarea=totalarea+sum(self.regions(i).Area); - end - else - totalarea=totalarea+sum(self.regions(region).Area); - end - - end % }}} - function [mpartition,npartition]=partition(self,varargin) % {{{ - mpartition=zeros(self.nregions(),1); - npartition=zeros(self.nregions(),1); + self=setdefaultparameters(self); + + options=pairoptions(varargin{:}); + + self.root=getfieldvalue(options,'root'); + region_names=getfieldvalue(options,'filenames'); + boxes_filename=getfieldvalue(options,'boxes'); + + %first read boxes regions shapefile: + disp('reading boxes for each region'); + self.boxes=shpread([self.root '/' boxes_filename]); + + %read the shape files and create the regions: counter=0; - for i=1:self.nregions(), - mpartition(i)=counter+1; - npartition(i)=counter+self.nglaciers(i); - counter=counter+self.nglaciers(i); - end - - end % }}} - 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(), - pos=find(gid>=mpartition(i) & gid<=npartition(i)); - rid(pos)=i; - lid(pos)=gid(pos)-mpartition(i)+1; - end - - end % }}} - function [gid]=lidtogid(self,rid,lid) % {{{ - [mpartition,npartition]=self.partition(); - gid=mpartition(rid)+lid-1; - end % }}} + self.regions=struct(); + for i=1:length(region_names), + disp(['reading region: ' region_names{i}]); + self.regions(i).name=region_names{i}; + self.regions(i).id=i; + contours=shpread([self.root '/' self.regions(i).name '.shp']); + + %we can't keep all the info: build arrays of centroids instead of reading + %the contours. + O1Region=zeros(length(contours),1); + O2Region=zeros(length(contours),1); + Area=zeros(length(contours),1); + CenLon=zeros(length(contours),1); + CenLat=zeros(length(contours),1); + for j=1:length(contours), + O1Region(j)=str2num(contours(j).O1Region); + O2Region(j)=str2num(contours(j).O2Region); + Area(j)=contours(j).Area; + CenLon(j)=contours(j).CenLon; + CenLat(j)=contours(j).CenLat; + end + self.regions(i).Area=Area; + self.regions(i).O1Region=O1Region; + self.regions(i).O2Region=O2Region; + self.regions(i).CenLat=CenLat; + self.regions(i).CenLon=CenLon; + self.regions(i).lids=[1:length(contours)]'; + self.regions(i).gids=self.regions(i).lids+counter; + counter=counter+length(contours); + end + end + %}}} function disp(self) % {{{ disp(sprintf(' Glacier inventory:')); @@ -162,75 +84,9 @@ end % }}} - function connectivity_plot(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); - - proj=laea(lat,long); - %plot([latel;latel(1)],[longel;longel(1)],'r-*') - %plot(lat,long,'k*'); - [latel;lat] - [longel;long] - [x,y]=gdaltransform([latel;lat],[longel;long],'EPSG:4326',proj) - error; - - end - - end % }}} - function checkconsistency(slm,solutiontype) % {{{ - - %is the coupler turned on? - for i=1:length(slm.icecaps), - if slm.icecaps{i}.transient.iscoupler==0, - warning(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name)); - end - end - - if slm.earth.transient.iscoupler==0, - warning('sealevelmodel checkconsistenty error: earth model should have the transient coupler option turned on!'); - end - - %check that the transition vectors have the right size: - for i=1:length(slm.icecaps), - if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.slr.transitions{i}), - error(['sealevelmodel checkconsistenty issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]); - end - end - - %check that run_frequency is the same everywhere: - for i=1:length(slm.icecaps), - if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency, - error(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name)); - end - end - - %make sure steric_rate is the same everywhere: - for i=1:length(slm.icecaps), - md= slm.icecaps{i}; - if ~isempty(find(md.slr.steric_rate - slm.earth.slr.steric_rate(slm.earth.slr.transitions{i}))), - error(sprintf('steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name)); - end - end - - end - %}}} function mesh_connectivity(self,mesh) % {{{ plotr=1; - %The way we run this is through the O2 zones defined in o2regions. We go through + %The way we run this is through the O2 zones defined in boxes. 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 @@ -263,6 +119,6 @@ %Go through O2 regions: - for i=1:length(self.o2regions), - string=self.o2regions(i).RGI_CODE; + for i=1:length(self.boxes), + string=self.boxes(i).RGI_CODE; disp(['progressing with region ' num2str(i) ' ' string]); offset=findstr(string,'-'); @@ -273,5 +129,5 @@ if ~isempty(glaciers), %find lat,long for laea projection: - box=self.o2regions(i).BoundingBox; + box=self.boxes(i).BoundingBox; long0=mean(box(:,1)); lat0=mean(box(:,2)); @@ -399,5 +255,81 @@ end % }}} - + function totalarea=area(self,varargin) % {{{ + region=-1; + totalarea=0; + if nargin==2, + region=varargin{1}; + end + if region==-1, + %figure out the areas of everybody: + for i=1:self.nregions(), + totalarea=totalarea+sum(self.regions(i).Area); + end + else + totalarea=totalarea+sum(self.regions(region).Area); + end + + end % }}} + function [mpartition,npartition]=partition(self,varargin) % {{{ + mpartition=zeros(self.nregions(),1); + npartition=zeros(self.nregions(),1); + counter=0; + for i=1:self.nregions(), + mpartition(i)=counter+1; + npartition(i)=counter+self.nglaciers(i); + counter=counter+self.nglaciers(i); + end + + end % }}} + function n = nregions(self) % {{{ + n=length(self.regions); + end + %}}} + function counter = nglaciers(self,varargin) % {{{ + + if nargin==1, + region=-1; %all regions. + else + region=varargin{1}; %only one. + end + + if region==-1, + counter=0; + for i=1:self.nregions(), + counter=counter+length(self.regions(i).Area); + end + else + counter=length(self.regions(region).Area); + end + end + %}}} + function name=ridtoname(self,rid) % {{{ + + fullname=self.region(rid).name; + name=fullname(10:end); + + end % }}} + function [gid]=lidtogid(self,rid,lid) % {{{ + [mpartition,npartition]=self.partition(); + gid=mpartition(rid)+lid-1; + end % }}} + 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(), + pos=find(gid>=mpartition(i) & gid<=npartition(i)); + rid(pos)=i; + lid(pos)=gid(pos)-mpartition(i)+1; + end + + end % }}} + function units(self) % {{{ + disp('Glacier inventory units: '); + disp(' areas: km^2'); + disp(' mass: Gt'); + + end + %}}} end end