Changeset 24156


Ignore:
Timestamp:
09/23/19 15:59:09 (5 years ago)
Author:
Eric.Larour
Message:

CHG: changes to work towards a better sea level model representation of bsins.

Location:
issm/trunk-jpl/src/m/classes
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/basin.m

    r24001 r24156  
    99                epsg              = 3426;
    1010                name              = '';
     11                continent         = '';
    1112        end
    1213        methods (Static)
    … …  
    2930                                        self.boundaries=getfieldvalue(options,'boundaries',{});
    3031                                        self.name=getfieldvalue(options,'name','');
     32                                        self.continent=getfieldvalue(options,'continent','');
    3133                                        self.epsg=getfieldvalue(options,'epsg',3426);
    3234                        end
    … …  
    3436                function self = setdefaultparameters(self) % {{{
    3537                        self.name='';
     38                        self.continent='';
    3639                        self.epsg=3426;
    3740                        self.boundaries={};
    … …  
    4043                function disp(self) % {{{
    4144                        disp(sprintf('   basin parameters:'));
     45                        fielddisplay(self,'continent','continent name');
    4246                        fielddisplay(self,'name','basin name');
    4347                        fielddisplay(self,'epsg','epsg projection number for the entire basin');
    … …  
    4751                        end
    4852
     53                end % }}}
     54                function boolean=isnameany(self,varargin) % {{{
     55                        boolean=0;
     56                        for  i=1:length(varargin),
     57                                if strcmpi(self.name,varargin{i}),
     58                                        boolean=1;
     59                                        break;
     60                                end
     61                        end
     62                end % }}}
     63                function boolean=iscontinentany(self,varargin) % {{{
     64                        boolean=0;
     65                        for  i=1:length(varargin),
     66                                if strcmpi(self.continent,varargin{i}),
     67                                        boolean=1;
     68                                        break;
     69                                end
     70                        end
    4971                end % }}}
    5072                function output=outputname(self,varargin) % {{{
    … …  
    87109                                boundary=self.boundaries{i};
    88110                                contour=boundary.edges();
    89                                 [contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.projection()),sprintf('EPSG:%i',self.epsg));
     111                                [contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.epsg),sprintf('EPSG:%i',self.epsg));
    90112                                x=[x;contour.x];
    91113                                y=[y;contour.y];
    92114                        end
     115                        %close the contour:
     116                        if x(end)~=x(1) | y(end)~=y(1),
     117                                x(end)=x(1); y(end)=y(1);
     118                        end
     119
    93120                        out.x=x;
    94121                        out.y=y;
     122                        out.nods=length(x);
    95123                end % }}}
    96124                function output=shapefilecrop(self,varargin) % {{{

  • issm/trunk-jpl/src/m/classes/boundary.m

    r24001 r24156  
    6666                end
    6767                end % }}}
    68                 function output=projection(self) % {{{
    69                         output=self.epsg;
    70                 end % }}}
    7168                function plot(self,varargin) % {{{
    7269                        %recover options

  • issm/trunk-jpl/src/m/classes/mesh3dsurface.m

    r22955 r24156  
    1717                long                        = NaN;
    1818                r                           = NaN;
     19                area                        = NaN;
    1920
    2021                vertexonboundary            = NaN;
    … …  
    115116                        fielddisplay(obj,'long','vertices longitude [degrees]');
    116117                        fielddisplay(obj,'r','vertices radius [m]');
     118                        fielddisplay(obj,'area','elemental areas [m^2]');
    117119                       
    118120                        fielddisplay(obj,'edges','edges of the 2d mesh (vertex1 vertex2 element1 element2)');
    … …  
    180182                        writejs1Darray(fid,[modelname '.mesh.long'],self.long);
    181183                        writejs1Darray(fid,[modelname '.mesh.r'],self.r);
     184                        writejs1Darray(fid,[modelname '.mesh.area'],self.area);
    182185                        writejs1Darray(fid,[modelname '.mesh.vertexonboundary'],self.vertexonboundary);
    183186                        writejs2Darray(fid,[modelname '.mesh.edges'],self.edges);

  • issm/trunk-jpl/src/m/classes/sealevelmodel.m

    r22955 r24156  
    99%      slm = sealevel('icecap',md_greenland,'icecap',md_antarctica,'earth',md_earth);
    1010
    11 classdef sealevelmodel
     11classdef sealevelmodel < handle
    1212        properties (SetAccess=public) %Model fields
    1313                % {{{
    1414                icecaps          = {}; % list of ice cap models
    1515                earth            = 0;  % model for the whole earth
     16                basins           = {}; % list  of basins, matching icecaps, where shapefile info is held.
    1617                cluster          = 0;
    1718                miscellaneous    = 0;
    … …  
    1920                private          = 0;
    2021                range            = 0;
    21                 mergedcaps     = 0;
     22                mergedcaps       = 0;
     23                transitions      = {};
    2224                %}}}
    2325        end
    … …  
    8789                        slm.cluster           = generic();
    8890                        slm.range             = {};
     91                        slm.transitions       = {};
    8992                end
    9093                %}}}
    … …  
    128131                                                        continue;
    129132                                                end
    130                                                                                         end
     133                                        end
    131134                                end
    132135                                self.mergedcaps{2*(i-1)+1}=md;
    … …  
    137140                                disp(sprintf('%i: %s',i,self.icecaps{i}.miscellaneous.name));
    138141                        end
     142                end % }}}
     143                function addbasin(self,bas) % {{{
     144                if ~strcmpi(class(bas),'basin')
     145                        error('addbasin method only takes a ''basin'' class object as input');
     146                end;
     147                self.basins{end+1}=bas;
     148                end % }}}
     149                function intersections(self,varargin) % {{{
     150
     151                options=pairoptions(varargin{:});
     152                force=getfieldvalue(options,'force',0);
     153
     154                for i=1:length(self.icecaps),
     155                        mdi=self.icecaps{i};
     156                        mdi=TwoDToThreeD(mdi);
     157               
     158                        disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
     159               
     160                        self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
     161
     162                end
     163
     164                end % }}}
     165                function checkintersections(self) % {{{
     166                flags=zeros(self.earth.mesh.numberofvertices,1);
     167                for i=1:length(self.basins),
     168                        flags(self.transitions{i})=i;
     169                end
     170                plotmodel(self.earth,'data',flags,'coastline','on');
     171
     172                end % }}}
     173                function baslist=basinindx(self,varargin) % {{{
     174                        options=pairoptions(varargin{:});
     175                        continent=getfieldvalue(options,'continent','all');
     176                        bas=getfieldvalue(options,'basin','all');
     177
     178                        %expand continent list: {{{
     179                        if iscell(continent),
     180                                if length(continent)==1,
     181                                         if strcmpi(continent{1},'all'),
     182                                                 %need to transform this into a list of continents:
     183                                                 continent={};
     184                                                 for i=1:length(self.basins),
     185                                                         continent{end+1}=self.basins{i}.continent;
     186                                                 end
     187                                                 continent=unique(continent);
     188                                         end
     189                                else
     190                                        %nothing to do, we have a list of continents
     191                                end
     192                        else
     193                                if strcmpi(continent,'all'),
     194                                        %need to transform this into a list of continents:
     195                                        continent={};
     196                                        for i=1:length(self.basins),
     197                                                 continent{end+1}=self.basins{i}.continent;
     198                                        end
     199                                        continent=unique(continent);
     200                                else
     201                                        continent={continent};
     202                                end
     203                        end
     204                        %}}}
     205                        %expand basins list using the continent list above and the extra bas discriminator: %{{{
     206                        if iscell(bas),
     207                                if length(bas)==1,
     208                                         if strcmpi(bas{1},'all'),
     209                                                 %need to transform this into a list of basins:
     210                                                 baslist=[];
     211                                                 for i=1:length(self.basins),
     212                                                         if self.basins{i}.iscontinentany(continent{:}),
     213                                                                 baslist(end+1)=i;
     214                                                         end
     215                                                 end
     216                                                 baslist=unique(baslist);
     217                                         else
     218                                                 bas=bas{1};
     219                                                 baslist=[];
     220                                                 for i=1:length(self.basins),
     221                                                         if self.basins{i}.iscontinentany(continent{:}),
     222                                                                 if self.basins{i}.isnameany(bas),
     223                                                                         baslist(end+1)=i;
     224                                                                 end
     225                                                         end
     226                                                 end
     227
     228                                         end
     229                                else
     230                                        %we have a list of basin names:
     231                                        baslist=[];
     232                                        for i=1:length(bas),
     233                                                basname=bas{i};
     234                                                for j=1:length(self.basins),
     235                                                        if self.basins{j}.iscontinentany(continent{:}),
     236                                                                if self.basins{j}.isnameany(basname),
     237                                                                        baslist(end+1)=j;
     238                                                                end
     239                                                        end
     240                                                end
     241                                                baslist=unique(baslist);
     242                                        end
     243                                end
     244                        else
     245                                if strcmpi(bas,'all'),
     246                                        baslist=[];
     247                                        for i=1:length(self.basins),
     248                                                if self.basins{i}.iscontinentany(continent{:}),
     249                                                        baslist(end+1)=i;
     250                                                end
     251                                        end
     252                                        baslist=unique(baslist);
     253                                else
     254                                        baslist=[];
     255                                        for i=1:length(self.basins),
     256                                                if self.basins{i}.iscontinentany(continent{:}),
     257                                                        if self.basins{i}.isnameany(bas),
     258                                                                baslist(end+1)=i;
     259                                                        end
     260                                                end
     261                                        end
     262                                        baslist=unique(baslist);
     263                                end
     264                        end
     265                        %}}}
     266
     267                end % }}}
     268                function addicecap(self,md) % {{{
     269                if ~strcmpi(class(md),'model')
     270                        error('addicecap method only takes a ''model'' class object as input');
     271                end
     272                self.icecaps{end+1}=md;
     273                end % }}}
     274                function basinsplot3d(self,varargin) % {{{
     275                for i=1:length(self.basins),
     276                        self.basins{i}.plot3d(varargin{:});
     277                end
     278                end % }}}
     279                function caticecaps(self,varargin) % {{{
     280                       
     281                        %recover options:
     282                        options=pairoptions(varargin{:});
     283                        tolerance=getfieldvalue(options,'tolerance',.65);
     284                        loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
     285       
     286                        %make 3D model:
     287                        models=self.icecaps;
     288                        for i=1:length(models),
     289                                models{i}=TwoDToThreeD(models{i});
     290                        end
     291                       
     292                        %Plug all models together:
     293                        md=models{1};
     294                        for i=2:length(models),
     295                                md=modelmerge3d(md,models{i},'tolerance',tolerance);
     296                                md.private.bamg.landmask=[md.private.bamg.landmask;models{i}.private.bamg.landmask];
     297                        end
     298
     299                        %Look for lone edges if asked for it: {{{
     300                        if loneedgesdetect,
     301                                edges=loneedges(md);
     302                                plotmodel(md,'data',md.mask.land_levelset);
     303                                hold on;
     304                                for i=1:length(edges),
     305                                        ind1=edges(i,1);
     306                                        ind2=edges(i,2);
     307                                        %plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'r*-');
     308                                        plot3([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],[md.mesh.z(ind1),md.mesh.z(ind2)],'g*-');
     309                                end
     310                        end %}}}
     311       
     312                        %Plug into earth:
     313                        self.earth=md;
     314
    139315                end % }}}
    140316                function viscousiterations(self) % {{{
Note: See TracChangeset for help on using the changeset viewer.