![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 24854
- Timestamp:
- 05/14/20 14:00:46 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
r24853 r24854 6 6 % where varargin is a variable list of options: 7 7 % 8 % Usage: 9 % rgi = glacier_inventory('root',shapefileroot,... 10 % 'filenames',region_names,... 11 % 'boxes', region_boxes); 8 12 % Example: 9 % rgi = glacier_inventory('shapefileroot',shapefileroot,'region_names',region_names); 13 % rgi = glacier_inventory('root','~/ModelData',... 14 % 'filenames',{'01_rgi60_Alaska','02_rgi60_WesternCanadaUS'},... 15 % 'boxes','00_rgi60_O2Regions.shp'); 10 16 % 11 % 17 % Watch out boxes do not have to match the region shapefiles. Example: O1regions vs O2 regions shapefiles. 12 18 13 19 classdef glacier_inventory < handle 14 20 properties (SetAccess=public) %Model fields 15 region_names = {}; 16 shapefileroot = ''; 21 root = ''; 17 22 regions = struct(); 18 o2regions= struct();23 boxes = struct(); 19 24 element_connectivity = []; 20 25 glacier_connectivity = []; … … 24 29 function self = glacier_inventory(varargin) % {{{ 25 30 26 if nargin==0, 27 self=setdefaultparameters(self); 28 else 29 self=setdefaultparameters(self); 30 31 options=pairoptions(varargin{:}); 32 33 self.shapefileroot=getfieldvalue(options,'shapefileroot'); 34 self.region_names=getfieldvalue(options,'region_names'); 35 o2regions_shapefile=getfieldvalue(options,'o2regions_shapefile'); 36 37 %first read o2 regions shapefile: 38 self.o2regions=shpread(o2regions_shapefile); 39 40 %read the shape files and create the regions: 41 counter=0; 42 self.regions=struct(); 43 for i=1:self.nregions(), 44 self.regions(i).name=self.region_names{i}; 45 disp(['reading region: ' self.regions(i).name]); 46 self.regions(i).id=i; 47 contours=shpread([self.shapefileroot '/' self.regions(i).name '.shp']); 48 49 %we can't keep all the info: build arrays of centroids: 50 O1Region=zeros(length(contours),1); 51 O2Region=zeros(length(contours),1); 52 Area=zeros(length(contours),1); 53 CenLon=zeros(length(contours),1); 54 CenLat=zeros(length(contours),1); 55 for j=1:length(contours), 56 O1Region(j)=str2num(contours(j).O1Region); 57 O2Region(j)=str2num(contours(j).O2Region); 58 Area(j)=contours(j).Area; 59 CenLon(j)=contours(j).CenLon; 60 CenLat(j)=contours(j).CenLat; 61 end 62 63 self.regions(i).Area=Area; 64 self.regions(i).O1Region=O1Region; 65 self.regions(i).O2Region=O2Region; 66 self.regions(i).CenLat=CenLat; 67 self.regions(i).CenLon=CenLon; 68 self.regions(i).lids=[1:length(contours)]'; 69 self.regions(i).gids=self.regions(i).lids+counter; 70 counter=counter+length(contours); 71 end 72 end 73 end 74 %}}} 75 function inv = setdefaultparameters(inv) % {{{ 76 end 77 %}}} 78 function n = nregions(self) % {{{ 79 n=length(self.region_names); 80 end 81 %}}} 82 function units(self) % {{{ 83 disp('Glacier inventory units: '); 84 disp(' areas: km^2'); 85 disp(' mass: Gt'); 86 87 end 88 %}}} 89 function counter = nglaciers(self,varargin) % {{{ 90 91 if nargin==1, 92 region=-1; %all regions. 93 else 94 region=varargin{1}; %only one. 95 end 96 97 if region==-1, 98 counter=0; 99 for i=1:self.nregions(), 100 counter=counter+length(self.regions(i).Area); 101 end 102 else 103 counter=length(self.regions(region).Area); 104 end 105 end 106 %}}} 107 function name=ridtoname(self,rid) % {{{ 108 109 fullname=self.region_names{rid}; 110 name=fullname(10:end); 111 112 end % }}} 113 function totalarea=area(self,varargin) % {{{ 114 region=-1; 115 totalarea=0; 116 if nargin==2, 117 region=varargin{1}; 118 end 119 if region==-1, 120 %figure out the areas of everybody: 121 for i=1:self.nregions(), 122 totalarea=totalarea+sum(self.regions(i).Area); 123 end 124 else 125 totalarea=totalarea+sum(self.regions(region).Area); 126 end 127 128 end % }}} 129 function [mpartition,npartition]=partition(self,varargin) % {{{ 130 mpartition=zeros(self.nregions(),1); 131 npartition=zeros(self.nregions(),1); 31 self=setdefaultparameters(self); 32 33 options=pairoptions(varargin{:}); 34 35 self.root=getfieldvalue(options,'root'); 36 region_names=getfieldvalue(options,'filenames'); 37 boxes_filename=getfieldvalue(options,'boxes'); 38 39 %first read boxes regions shapefile: 40 disp('reading boxes for each region'); 41 self.boxes=shpread([self.root '/' boxes_filename]); 42 43 %read the shape files and create the regions: 132 44 counter=0; 133 for i=1:self.nregions(), 134 mpartition(i)=counter+1; 135 npartition(i)=counter+self.nglaciers(i); 136 counter=counter+self.nglaciers(i); 137 end 138 139 end % }}} 140 function [lid,rid]=gidtolid(self,gid) % {{{ 141 [mpartition,npartition]=self.partition(); 142 lid=zeros(length(gid),1); 143 rid=zeros(length(gid),1); 144 for i=1:self.nregions(), 145 pos=find(gid>=mpartition(i) & gid<=npartition(i)); 146 rid(pos)=i; 147 lid(pos)=gid(pos)-mpartition(i)+1; 148 end 149 150 end % }}} 151 function [gid]=lidtogid(self,rid,lid) % {{{ 152 [mpartition,npartition]=self.partition(); 153 gid=mpartition(rid)+lid-1; 154 end % }}} 45 self.regions=struct(); 46 for i=1:length(region_names), 47 disp(['reading region: ' region_names{i}]); 48 self.regions(i).name=region_names{i}; 49 self.regions(i).id=i; 50 contours=shpread([self.root '/' self.regions(i).name '.shp']); 51 52 %we can't keep all the info: build arrays of centroids instead of reading 53 %the contours. 54 O1Region=zeros(length(contours),1); 55 O2Region=zeros(length(contours),1); 56 Area=zeros(length(contours),1); 57 CenLon=zeros(length(contours),1); 58 CenLat=zeros(length(contours),1); 59 for j=1:length(contours), 60 O1Region(j)=str2num(contours(j).O1Region); 61 O2Region(j)=str2num(contours(j).O2Region); 62 Area(j)=contours(j).Area; 63 CenLon(j)=contours(j).CenLon; 64 CenLat(j)=contours(j).CenLat; 65 end 66 self.regions(i).Area=Area; 67 self.regions(i).O1Region=O1Region; 68 self.regions(i).O2Region=O2Region; 69 self.regions(i).CenLat=CenLat; 70 self.regions(i).CenLon=CenLon; 71 self.regions(i).lids=[1:length(contours)]'; 72 self.regions(i).gids=self.regions(i).lids+counter; 73 counter=counter+length(contours); 74 end 75 end 76 %}}} 155 77 function disp(self) % {{{ 156 78 disp(sprintf(' Glacier inventory:')); … … 162 84 163 85 end % }}} 164 function connectivity_plot(self,md) % {{{165 disp(sprintf(' Checking glacier connectivity:'));166 figure(1),clf,hold on;167 168 %some quick checks: % {{{169 if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)),170 error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!');171 end172 if ~isempty(find(self.glacier_connectivity<0)),173 error('glacier_connectivity table is pointing to elements that are < 0!');174 end %}}}175 176 for gid=1:length(self.glacier_connectivity),177 el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:));178 [lid,rid]=self.gidtolid(gid);179 lat=self.regions(rid).CenLat(lid);180 long=self.regions(rid).CenLon(lid);181 182 proj=laea(lat,long);183 %plot([latel;latel(1)],[longel;longel(1)],'r-*')184 %plot(lat,long,'k*');185 [latel;lat]186 [longel;long]187 [x,y]=gdaltransform([latel;lat],[longel;long],'EPSG:4326',proj)188 error;189 190 end191 192 end % }}}193 function checkconsistency(slm,solutiontype) % {{{194 195 %is the coupler turned on?196 for i=1:length(slm.icecaps),197 if slm.icecaps{i}.transient.iscoupler==0,198 warning(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));199 end200 end201 202 if slm.earth.transient.iscoupler==0,203 warning('sealevelmodel checkconsistenty error: earth model should have the transient coupler option turned on!');204 end205 206 %check that the transition vectors have the right size:207 for i=1:length(slm.icecaps),208 if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.slr.transitions{i}),209 error(['sealevelmodel checkconsistenty issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]);210 end211 end212 213 %check that run_frequency is the same everywhere:214 for i=1:length(slm.icecaps),215 if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency,216 error(sprintf('sealevelmodel checkconsistenty error: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));217 end218 end219 220 %make sure steric_rate is the same everywhere:221 for i=1:length(slm.icecaps),222 md= slm.icecaps{i};223 if ~isempty(find(md.slr.steric_rate - slm.earth.slr.steric_rate(slm.earth.slr.transitions{i}))),224 error(sprintf('steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));225 end226 end227 228 end229 %}}}230 86 function mesh_connectivity(self,mesh) % {{{ 231 87 232 88 plotr=1; 233 89 234 %The way we run this is through the O2 zones defined in o2regions. We go through90 %The way we run this is through the O2 zones defined in boxes. We go through 235 91 %these regions, figure out the centroid, figure out how many elements are close to 236 92 %this centroid (very important to do this through vertex lat,long, and not through elemnet … … 263 119 264 120 %Go through O2 regions: 265 for i=1:length(self. o2regions),266 string=self. o2regions(i).RGI_CODE;121 for i=1:length(self.boxes), 122 string=self.boxes(i).RGI_CODE; 267 123 disp(['progressing with region ' num2str(i) ' ' string]); 268 124 offset=findstr(string,'-'); … … 273 129 if ~isempty(glaciers), 274 130 %find lat,long for laea projection: 275 box=self. o2regions(i).BoundingBox;131 box=self.boxes(i).BoundingBox; 276 132 long0=mean(box(:,1)); 277 133 lat0=mean(box(:,2)); … … 399 255 400 256 end % }}} 401 257 function totalarea=area(self,varargin) % {{{ 258 region=-1; 259 totalarea=0; 260 if nargin==2, 261 region=varargin{1}; 262 end 263 if region==-1, 264 %figure out the areas of everybody: 265 for i=1:self.nregions(), 266 totalarea=totalarea+sum(self.regions(i).Area); 267 end 268 else 269 totalarea=totalarea+sum(self.regions(region).Area); 270 end 271 272 end % }}} 273 function [mpartition,npartition]=partition(self,varargin) % {{{ 274 mpartition=zeros(self.nregions(),1); 275 npartition=zeros(self.nregions(),1); 276 counter=0; 277 for i=1:self.nregions(), 278 mpartition(i)=counter+1; 279 npartition(i)=counter+self.nglaciers(i); 280 counter=counter+self.nglaciers(i); 281 end 282 283 end % }}} 284 function n = nregions(self) % {{{ 285 n=length(self.regions); 286 end 287 %}}} 288 function counter = nglaciers(self,varargin) % {{{ 289 290 if nargin==1, 291 region=-1; %all regions. 292 else 293 region=varargin{1}; %only one. 294 end 295 296 if region==-1, 297 counter=0; 298 for i=1:self.nregions(), 299 counter=counter+length(self.regions(i).Area); 300 end 301 else 302 counter=length(self.regions(region).Area); 303 end 304 end 305 %}}} 306 function name=ridtoname(self,rid) % {{{ 307 308 fullname=self.region(rid).name; 309 name=fullname(10:end); 310 311 end % }}} 312 function [gid]=lidtogid(self,rid,lid) % {{{ 313 [mpartition,npartition]=self.partition(); 314 gid=mpartition(rid)+lid-1; 315 end % }}} 316 function [lid,rid]=gidtolid(self,gid) % {{{ 317 [mpartition,npartition]=self.partition(); 318 lid=zeros(length(gid),1); 319 rid=zeros(length(gid),1); 320 for i=1:self.nregions(), 321 pos=find(gid>=mpartition(i) & gid<=npartition(i)); 322 rid(pos)=i; 323 lid(pos)=gid(pos)-mpartition(i)+1; 324 end 325 326 end % }}} 327 function units(self) % {{{ 328 disp('Glacier inventory units: '); 329 disp(' areas: km^2'); 330 disp(' mass: Gt'); 331 332 end 333 %}}} 402 334 end 403 335 end
Note:
See TracChangeset
for help on using the changeset viewer.
