![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 24844
- Timestamp:
- 05/11/20 19:16:40 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
r24843 r24844 16 16 shapefileroot = ''; 17 17 regions = struct(); 18 element_connectivity = []; 19 glacier_connectivity = []; 18 20 end 19 21 methods … … 73 75 end 74 76 %}}} 77 function counter = nglaciers(self,varargin) % {{{ 78 79 if nargin==1, 80 region=-1; %all regions. 81 else 82 region=varargin{1}; %only one. 83 end 84 85 if region==-1, 86 counter=0; 87 for i=1:self.nregions(), 88 counter=counter+length(self.regions(i).Area); 89 end 90 else 91 counter=length(self.regions(region).Area); 92 end 93 end 94 %}}} 95 function [mpartition,npartition]=partition(self,varargin) % {{{ 96 mpartition=zeros(self.nregions(),1); 97 npartition=zeros(self.nregions(),1); 98 counter=0; 99 for i=1:self.nregions(), 100 mpartition(i)=counter+1; 101 npartition(i)=counter+self.nglaciers(i); 102 counter=counter+self.nglaciers(i); 103 end 104 105 end % }}} 106 function [lid,rid]=gidtolid(self,gid) % {{{ 107 [mpartition,npartition]=self.partition(); 108 for i=1:self.nregions(), 109 if (gid>=mpartition(i) && gid <=npartition(i)), 110 rid=i; 111 lid=gid-mpartition(i)+1; 112 break; 113 end 114 end 115 116 end % }}} 117 function [gid]=lidtogid(self,rid,lid) % {{{ 118 [mpartition,npartition]=self.partition(); 119 gid=mpartition(rid)+lid-1; 120 end % }}} 75 121 function disp(self) % {{{ 76 122 disp(sprintf(' Glacier inventory:')); … … 79 125 for i=1:self.nregions(), 80 126 disp(sprintf(' region %i: ''%s'' %i glaciers ',i,self.regions(i).name,length(self.regions(i).Area))); 127 end 128 129 end % }}} 130 function connectivity_check(self,md) % {{{ 131 disp(sprintf(' Checking glacier connectivity:')); 132 figure(1),clf,hold on; 133 counter=0; 134 135 %some quick checks: % {{{ 136 if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)), 137 error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!'); 138 end 139 if ~isempty(find(self.glacier_connectivity<0)), 140 error('glacier_connectivity table is pointing to elements that are < 0!'); 141 end %}}} 142 143 for gid=1:length(self.glacier_connectivity), 144 el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:)); 145 x1=latel(1); x2=latel(2); x3=latel(3); 146 y1=longel(1); y2=longel(2); y3=longel(3); 147 148 [lid,rid]=self.gidtolid(gid); 149 x0=self.regions(rid).CenLat(lid); 150 y0=self.regions(rid).CenLon(lid); 151 152 if ~isintriangle(x0,x1,x2,x3,y0,y1,y2,y3), 153 %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*'); 154 %error(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el)); 155 disp(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el)); 156 counter=counter+1 157 end 158 end 159 160 end % }}} 161 function connectivity_check2(self,md) % {{{ 162 disp(sprintf(' Checking glacier connectivity:')); 163 figure(1),clf,hold on; 164 165 %some quick checks: % {{{ 166 if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)), 167 error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!'); 168 end 169 if ~isempty(find(self.glacier_connectivity<0)), 170 error('glacier_connectivity table is pointing to elements that are < 0!'); 171 end %}}} 172 173 for gid=1:length(self.glacier_connectivity), 174 el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:)); 175 [lid,rid]=self.gidtolid(gid); 176 lat=self.regions(rid).CenLat(lid); 177 long=self.regions(rid).CenLon(lid); 178 179 proj=laea(lat,long); 180 [longel;long] 181 [x,y]=gdaltransform([longel;long],[latel;lat],'EPSG:4326',proj); 182 x1=x(1); y1=y(1); x2=x(2); y2=y(2); x3=x(3); y3=y(3); x0=x(4); y0=y(4); 183 184 %test if x0,y0 is in the 1,2,3 triangle: 185 epsilon=10^-10; 186 delta=x2*y3-y2*x3-x1*y3+y1*x3+x1*y2-y1*x2; 187 188 %first area coordinate 189 area_1=((y3-y2)*(x0-x3)-(x3-x2)*(y0-y3))./delta 190 %second area coordinate 191 area_2=((x3-x1)*(y0-y3)-(y3-y1)*(x0-x3))./delta 192 %third area coordinate 193 area_3=1-area_1-area_2 194 195 isintriangle=((area_1>=0-epsilon) & (area_2>=0-epsilon) & (area_3>=0-epsilon)); 196 197 plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-'); 198 plot(x0,y0,'r*'); 199 if ~isintriangle, 200 error(); 201 end 202 203 end 204 205 end % }}} 206 function connectivity_plot(self,md) % {{{ 207 disp(sprintf(' Checking glacier connectivity:')); 208 figure(1),clf,hold on; 209 210 %some quick checks: % {{{ 211 if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)), 212 error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!'); 213 end 214 if ~isempty(find(self.glacier_connectivity<0)), 215 error('glacier_connectivity table is pointing to elements that are < 0!'); 216 end %}}} 217 218 for gid=1:length(self.glacier_connectivity), 219 el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:)); 220 [lid,rid]=self.gidtolid(gid); 221 lat=self.regions(rid).CenLat(lid); 222 long=self.regions(rid).CenLon(lid); 223 224 proj=laea(lat,long); 225 %plot([latel;latel(1)],[longel;longel(1)],'r-*') 226 %plot(lat,long,'k*'); 227 [latel;lat] 228 [longel;long] 229 [x,y]=gdaltransform([latel;lat],[longel;long],'EPSG:4326',proj) 230 error; 231 81 232 end 82 233 … … 119 270 end 120 271 %}}} 272 function mesh_connectivity(self,mesh) % {{{ 273 274 %initialize glacier_connectivity: 275 glacier_connectivity=zeros(self.nglaciers,1); 276 277 %assume maximum width for connectivity table and initialize 278 %as sparse matrix: 279 ny=3000; 280 element_connectivity=sparse(mesh.numberofelements,ny); 281 282 disp('Building element and glacier connectivity table: '); 283 284 partition_counter=0; 285 for i=1:self.nregions(), 286 region=self.regions(i); 287 288 disp([' region ' num2str(i) ': ' region.name]); 289 290 %reproject on a lambert centered on the region: 291 latel=mesh.lat(mesh.elements)*[1;1;1]/3; 292 longel=mesh.long(mesh.elements)*[1;1;1]/3; 293 for j=1:length(region.CenLat), 294 if mod(j,100)==0, 295 disp(num2str(j/length(region.CenLat)*100)); 296 end 297 298 lat0=region.CenLat(j); long0=region.CenLon(j); 299 %proj=laea(lat0,long0); 300 301 %to avoid folding of projections from behind the Earth, restrict ourselves 302 %to only locations around a lat,long radius around lat0,long0 303 list_elements=flaglatlongradius(latel,longel,lat0,long0,5); 304 305 %go through these elements and figure out if the glacier is inside: 306 found=0; 307 for k=1:length(list_elements), 308 el=list_elements(k); 309 lat=mesh.lat(mesh.elements(el,:)); 310 long=mesh.long(mesh.elements(el,:)); 311 x0=lat0; x1=lat(1); x2=lat(2); x3=lat(3); 312 y0=long0; y1=long(1); y2=long(2); y3=long(3); 313 314 if isintriangle(x0,x1,x2,x3,y0,y1,y2,y3), 315 %we have found the element where this glacier belongs: 316 found=1; 317 break; 318 end 319 end; 320 if ~found, %try again with a projection and all elements: % {{{ 321 proj=laea(lat0,long0); 322 list_elements=(1:mesh.numberofelements)'; 323 found=0; 324 [xp,yp]=gdaltransform([mesh.long;long0],[mesh.lat;lat0],'EPSG:4326',proj); 325 x0=xp(end); y0=yp(end); 326 for k=1:length(list_elements), 327 el=list_elements(k); 328 x=xp(mesh.elements(el,:)); y=yp(mesh.elements(el,:)); 329 x1=x(1); x2=x(2); x3=x(3); 330 y1=y(1); y2=y(2); y3=y(3); 331 if isintriangle(x0,x1,x2,x3,y0,y1,y2,y3), 332 %we have found the element where this glacier belongs: 333 found=1; 334 break; 335 end 336 end; 337 end % }}} 338 if ~found, %not found, plot %{{{ 339 figure(1),clf,hold on; 340 plot(x0,y0,'r*'); 341 for k=1:length(list_elements), 342 el=list_elements(k); 343 x=xp(mesh.elements(el,:)); y=yp(mesh.elements(el,:)); 344 x1=x(1); x2=x(2); x3=x(3); 345 y1=y(1); y2=y(2); y3=y(3); 346 plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-'); 347 end 348 error('cound not find an element, believe it or not'); %}}} 349 end 350 glacier_connectivity(partition_counter+j)=el; 351 end 352 partition_counter=partition_counter+length(region.CenLat); 353 end 354 355 %Now that we have all the elements for each glacier, we can build a 356 %connectivity table that is the reverse of the glacier connectivity table. 357 %For each element, this will tell us which glaciers belong to this element. 358 for j=1:length(glacier_connectivity), 359 el=glacier_connectivity(j); 360 if el==0, 361 continue; 362 end 363 count=element_connectivity(el,ny); 364 if count>ny, 365 error('need to enlarge connectivity table'); 366 end 367 element_connectivity(el,count+1)=partition_counter+j; 368 element_connectivity(el,ny)=count+1; 369 end 370 371 %Reduce the number of columns (we did not initially, so we chose an arbitrary ny: 372 nmax=max(element_connectivity(:,end)); 373 element_connectivity=element_connectivity(:,[1:nmax,ny]); 374 375 %Assign pointers: 376 self.element_connectivity=element_connectivity; 377 self.glacier_connectivity=glacier_connectivity; 378 379 end % }}} 121 380 end 122 381 end
Note:
See TracChangeset
for help on using the changeset viewer.
