Changeset 24779


Ignore:
Timestamp:
05/04/20 10:03:33 (5 years ago)
Author:
Eric.Larour
Message:

CHG: boundary and basin now rely on a proj string instead of an epsg code. Adds
flexibility (when combined for example with laea.m routine). Transported implied
changes into sealevelmodel, TwoDToThreeD and meshintersect3d.

sealevelmodel is also now depending on an extrat 'planet' field, so that
we can recover the right radius in TwoDToThreeD. This will pave the way for
planetary bodies application such as 'Europa', 'Titan', etc...

Location:
issm/trunk-jpl/src/m
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r24761 r24779  
    77        properties (SetAccess=public)
    88                boundaries        = {};
    9                 epsg              = 3426;
    109                name              = '';
    1110                continent         = '';
     11                proj              = epsg2proj(4326);
    1212        end
    1313        methods (Static)
    … …  
    3131                                        self.name=getfieldvalue(options,'name','');
    3232                                        self.continent=getfieldvalue(options,'continent','');
    33                                         self.epsg=getfieldvalue(options,'epsg',3426);
     33
     34                                        %figure out projection string:
     35                                        if exist(options,'epsg'),
     36                                                if exist(options,'proj'),
     37                                                        error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
     38                                                end
     39                                                epsg=getfieldvalue(options,'epsg');
     40                                                %convert into proj4:
     41                                                proj=epsg2proj(epsg);
     42                                        elseif exist(options,'proj'),
     43                                                if exist(options,'epsg'),
     44                                                        error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
     45                                                end
     46                                                proj=getfieldvalue(options,'proj');
     47                                        else
     48                                                proj= '+proj=longlat +datum=WGS84 +no_defs';
     49                                        end
     50                                        self.proj=proj;
    3451                        end
    3552                end % }}}
    … …  
    3754                        self.name='';
    3855                        self.continent='';
    39                         self.epsg=3426;
     56                        self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
    4057                        self.boundaries={};
    4158
    … …  
    4562                        fielddisplay(self,'continent','continent name');
    4663                        fielddisplay(self,'name','basin name');
    47                         fielddisplay(self,'epsg','epsg projection number for the entire basin');
     64                        fielddisplay(self,'proj','basin projection string (follows PROJ.4 standard)');
    4865                        fielddisplay(self,'boundaries','list of boundary objects');
    4966                        for i=1:length(self.boundaries),
    … …  
    87104                        %add option:
    88105                        for i=1:length(self.boundaries),
    89                                 self.boundaries{i}.plot('epsg',self.epsg,varargin{:});
     106                                self.boundaries{i}.plot('proj',self.proj,varargin{:});
    90107                        end
    91108
    … …  
    109126                                boundary=self.boundaries{i};
    110127                                contour=boundary.edges();
    111                                 [contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.epsg),sprintf('EPSG:%i',self.epsg));
     128                                [contour.x,contour.y]=gdaltransform(contour.x,contour.y,boundary.proj,self.proj);
    112129                                x=[x;contour.x];
    113130                                y=[y;contour.y];
    … …  
    122139                        out.nods=length(x);
    123140                end % }}}
     141                function checkconsistency(self,varargin) % {{{
     142               
     143                        %recover options
     144                        options=pairoptions(varargin{:});
     145
     146                        %figure out if two boundaries are identical:
     147                        for i=1:length(self.boundaries),
     148                                namei=self.boundaries{i}.shpfilename;
     149                                for j=i+1:length(self.boundaries),
     150                                        namej=self.boundaries{j}.shpfilename;
     151                                        if strcmpi(namei,namej),
     152                                                error(['Basin ' self.name ' has two identical boundaries named ' namei ]);
     153                                        end
     154                                end
     155                        end
     156
     157                        %go through boundaries and check their consistency:
     158                        for i=1:length(self.boundaries),
     159                                boundary=self.boundaries{i};
     160                                boundary.checkconsistency();
     161                        end
     162                end % }}}
    124163                function contourplot(self,varargin) % {{{
    125164                        contour=self.contour();
    … …  
    134173                        inshapefile=getfieldvalue(options,'shapefile');
    135174                        outputshapefile=getfieldvalue(options,'outputshapefile','');
    136                         epsgshapefile=getfieldvalue(options,'epsgshapefile');
     175
     176                        if (exist(options,'epsgshapefile') &  exist(options,'projshapefile')),
     177                                error('Basin shapefilecrop error message: cannot specify epsgshapefile and projshapefile at the same time');
     178                        end
     179                        if exist(options,'epsgshapefile'),
     180                                projshapefile=epsg2proj(getfieldvalue(options,'epsgshapefile'));
     181                        elseif exist(options,'projshapefile'),
     182                                projshapefile=getfieldvalue(options,'projshapefile');
     183                        else
     184                                error('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified');
     185                        end
     186
    137187
    138188                        %create list of contours that have critical length > threshold:  (in lat,long)
    … …  
    150200
    151201                        %project onto reference frame:
    152                         if self.epsg~=epsgshapefile,
    153                                 for i=1:length(contours),
    154                                         h=contours(i);
    155                                         [h.x,h.y]=gdaltransform(h.x,h.y,sprintf('EPSG:%i',epsgshapefile),sprintf('EPSG:%i',self.epsg));
    156                                         contours(i).x=h.x;
    157                                         contours(i).y=h.y;
    158                                 end
     202                        for i=1:length(contours),
     203                                h=contours(i);
     204                                [h.x,h.y]=gdaltransform(h.x,h.y,projshapefile,self.proj);
     205                                contours(i).x=h.x;
     206                                contours(i).y=h.y;
    159207                        end
    160208

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

    r24765 r24779  
    99                shpfilename       = '';
    1010                orientation       = 'normal';  %other value is 'reverse'
    11                 epsg              = 4326; %EPSG number, default value is WGS 84 Lat,Long reference frame.
     11                proj              = epsg2proj(4326); %Proj4.0  projection string. , default is WGS 84 (corresponds to epsg 4326)
     12
    1213        end
    1314        methods (Static)
    … …  
    3031                                        self.shpfilename=getfieldvalue(options,'shpfilename','');
    3132                                        self.orientation=getfieldvalue(options,'orientation','normal');
    32                                         self.epsg=getfieldvalue(options,'epsg',4326);
     33
     34                                        %figure out projection string:
     35                                        if exist(options,'epsg'),
     36                                                if exist(options,'proj'),
     37                                                        error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
     38                                                end
     39                                                epsg=getfieldvalue(options,'epsg');
     40                                                %convert into proj4:
     41                                                proj=epsg2proj(epsg);
     42                                        elseif exist(options,'proj'),
     43                                                if exist(options,'epsg'),
     44                                                        error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
     45                                                end
     46                                                proj=getfieldvalue(options,'proj');
     47                                        else
     48                                                proj= '+proj=longlat +datum=WGS84 +no_defs';
     49                                        end
     50                                        self.proj=proj;
    3351                        end
    3452                end % }}}
    … …  
    3755                self.shpfilename='';
    3856                self.orientation='normal';
    39                 self.epsg=4326;
     57                self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
    4058                end % }}}
    4159                function disp(self) % {{{
    … …  
    4563                        fielddisplay(self,'shpfilename','shape file name');
    4664                        fielddisplay(self,'orientation','orientation (default is ''normal'', can be ''reverse'')');
    47                         fielddisplay(self,'epsg','EPSG number defining projection for the shape file');
     65                        fielddisplay(self,'proj','shape file projection string (follows PROJ.4 standard)');
    4866
    4967                end % }}}
    … …  
    7795                        markersize=getfieldvalue(options,'markersize',1);
    7896                        unitmultiplier=getfieldvalue(options,'unit',1);
    79                         epsg=getfieldvalue(options,'epsg',4326);
    8097                        radius=getfieldvalue(options,'radius',6371012);
    8198                        aboveground=getfieldvalue(options,'aboveground',1000);
    … …  
    84101                        label=getfieldvalue(options,'label','on');
    85102
     103                        if (exist(options,'epsg') &  exist(options,'proj')),
     104                                error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
     105                        end
     106                        if exist(options,'epsg'),
     107                                proj=epsg2proj(getfieldvalue(options,'epsg'));
     108                        elseif exist(options,'proj'),
     109                                proj=getfieldvalue(options,'proj');
     110                        else
     111                                proj=epsg2proj(4326);
     112                        end
     113
    86114                        %read domain:
    87115                        [path,name,ext]=fileparts(self.shpfilename);
    … …  
    92120
    93121                        %convert boundary to another reference frame:  {{{
    94 
    95122                        for i=1:length(domain),
    96123                                try,
    97                                         [x,y] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),sprintf('EPSG:%i',epsg));
     124                                        [x,y] = gdaltransform(domain(i).x,domain(i).y,self.proj,proj);
    98125                                catch me
    99126                                        disp(me.message());
    … …  
    125152                        %}}}
    126153                end % }}}
     154                function checkconsistency(self,varargin) % {{{
     155                        %recover options
     156                        options=pairoptions(varargin{:});
     157                        tolerance=getfieldvalue(options,'tolerance',1e-5);
     158               
     159                        %recover edges:
     160                        edges=self.edges();
     161
     162                        if strcmpi(edges.Geometry,'Point'),
     163                                return;
     164                        else
     165                                %check we don't have two identical vertices:
     166                                x=edges.x; y=edges.y;
     167                                distance=sqrt( (x(1:end-1)-x(2:end)).^2+(y(1:end-1)-y(2:end)).^2);
     168                                dmax=max(distance);
     169                                isidentical=find(distance<dmax*tolerance);
     170                                if ~isempty(isidentical),
     171                                        error(['boundary ' shpfilename ' has two  vertices extremely close to one another']);
     172                                end
     173                        end
     174                end % }}}
    127175                function plot3d(self,varargin) % {{{
    128176                        %recover options
    … …  
    136184                        markersize=getfieldvalue(options,'markersize',1);
    137185                        unitmultiplier=getfieldvalue(options,'unit',1);
    138                         epsg=getfieldvalue(options,'epsg',4326);
    139186                        radius=getfieldvalue(options,'radius',6371012);
    140187                        aboveground=getfieldvalue(options,'aboveground',1000);
    … …  
    150197                        domain=shpread([self.shppath '/' name '.' ext]);
    151198
    152                         if epsg==4326,
    153                                 %convert boundary to lat,long: {{{
    154 
    155                                 for i=1:length(domain),
    156                                         try,
    157                                                 [lat,long] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),'EPSG:4326');
    158                                         catch me
    159                                                 disp(me.message());
    160                                                 self.disp();
    161                                         end
    162                                         domain(i).x=long; domain(i).y=lat;
    163                                 end
    164 
    165                                 for i=1:length(domain),
    166 
    167                                         %make sure lat,long are what they are supposed to be:
    168                                         %if any(domain(i).x>90 | domain(i).x<-90),
    169                                         %       long=domain(i).x; lat=domain(i).y;
    170                                         %else
    171                                         %       long=domain(i).y; lat=domain(i).x;
    172                                         %end
    173 
    174                                         %project on x,y,z reference frame.
    175                                         [x,y,z]=AboveGround(domain(i).x,domain(i).y,radius,aboveground);
    176                                         [xt,yt,zt]=AboveGround(domain(i).x+offset,domain(i).y+offset,radius,300000);
    177                                         hold on;
    178                                         if length(x)==1,
    179                                                 p=plot3(x,y,z,'k*');
    180                                                 set(p,'MarkerSize',markersize);
    181                                                 if strcmpi(label,'on'),
    182                                                         t=text(xt,yt,zt,self.shpfilename,'FontSize',fontsize);
    183                                                 end
    184                                         else
    185                                                 p=plot3(x,y,z,'k-');
    186                                                 mid=floor(length(xt)/2);
    187                                                 if strcmpi(label,'on'),
    188                                                         text(xt(mid),yt(mid),zt(mid),self.shpfilename,'FontSize',fontsize);
    189                                                 end
    190                                         end
    191                                         set(p,'Color',color);
    192                                         set(p,'LineWidth',linewidth);
    193 
    194                                 end
    195                                 %}}}
    196                         else
    197                                 %convert boundary to another reference frame:  {{{
    198 
    199                                 for i=1:length(domain),
    200                                         try,
    201                                                 [x,y] = gdaltransform(domain(i).x,domain(i).y,sprintf('EPSG:%i',self.epsg),sprintf('EPSG:%i',epsg));
    202                                         catch me
    203                                                 disp(me.message());
    204                                                 self.disp();
    205                                         end
    206                                         domain(i).x=x; domain(i).y=y;
    207                                 end
    208 
    209                                 for i=1:length(domain),
    210                                         hold on;
    211                                         if length(x)==1,
    212                                                 p=plot(x,y,'k*');
    213                                                 set(p,'MarkerSize',markersize);
    214                                                 if strcmpi(label,'on'),
    215                                                         t=text(x,y,self.shpfilename,'FontSize',fontsize);
    216                                                 end
    217                                         else
    218                                                 p=plot(x,y,'k-');
    219                                                 if strcmpi(label,'on'),
    220                                                         text(sum(x)/length(x),sum(y)/length(y),self.shpfilename,'FontSize',fontsize);
    221                                                 end
    222                                         end
    223                                         set(p,'Color',color);
    224                                         set(p,'LineWidth',linewidth);
    225                                 end
    226                                 %}}}
    227                         end
     199                        for i=1:length(domain),
     200                                try,
     201                                        [lat,long] = gdaltransform(domain(i).x,domain(i).y,self.proj,'EPSG:4326');
     202                                catch me
     203                                        disp(me.message());
     204                                        self.disp();
     205                                end
     206                                domain(i).x=long; domain(i).y=lat;
     207                        end
     208
     209                        for i=1:length(domain),
     210
     211                                %project on x,y,z reference frame.
     212                                [x,y,z]=AboveGround(domain(i).x,domain(i).y,radius,aboveground);
     213                                [xt,yt,zt]=AboveGround(domain(i).x+offset,domain(i).y+offset,radius,300000);
     214                                hold on;
     215                                if length(x)==1,
     216                                        p=plot3(x,y,z,'k*');
     217                                        set(p,'MarkerSize',markersize);
     218                                        if strcmpi(label,'on'),
     219                                                t=text(xt,yt,zt,self.shpfilename,'FontSize',fontsize);
     220                                        end
     221                                else
     222                                        p=plot3(x,y,z,'k-');
     223                                        mid=floor(length(xt)/2);
     224                                        if strcmpi(label,'on'),
     225                                                text(xt(mid),yt(mid),zt(mid),self.shpfilename,'FontSize',fontsize);
     226                                        end
     227                                end
     228                                set(p,'Color',color);
     229                                set(p,'LineWidth',linewidth);
     230
     231                        end
     232
    228233                end % }}}
    229234        end

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

    r24728 r24779  
    2323                transitions      = {};
    2424                eltransitions      = {};
     25                planet           = '';
    2526                %}}}
    2627        end
    … …  
    3940                               
    4041                                %recover the earth model:
    41                                 slm.earth = getfieldvalue(options,'earth');
     42                                slm.earth = getfieldvalue(options,'earth',0);
     43
     44                                %set planet type:
     45                                slm.planet=getfieldvalue(options,'planet','earth');
     46
    4247                        end
    4348                end
    … …  
    9196                        slm.range             = {};
    9297                        slm.transitions       = {};
    93                         slm.eltransitions       = {};
     98                        slm.eltransitions     = {};
     99                        slm.planet            = 'earth';
    94100                end
    95101                %}}}
    … …  
    153159                options=pairoptions(varargin{:});
    154160                force=getfieldvalue(options,'force',0);
     161               
     162                %initialize, to avoid issues of having more transitions than meshes.
     163                self.transitions={}; self.eltransitions={};
    155164
    156165                %for elements:
    … …  
    161170                for i=1:length(self.icecaps),
    162171                        mdi=self.icecaps{i};
    163                         mdi=TwoDToThreeD(mdi);
     172                        mdi=TwoDToThreeD(mdi,self.planet);
    164173       
    165174                        %for elements:
    … …  
    182191                end
    183192                plotmodel(self.earth,'data',flags,'coastline','on');
     193
     194                end % }}}
     195                function checkbasinconsistency(self) % {{{
     196                        for i=1:length(self.basins),
     197                                self.basins{i}.checkconsistency();
     198                        end
    184199
    185200                end % }}}
    … …  
    300315                        models=self.icecaps;
    301316                        for i=1:length(models),
    302                                 models{i}=TwoDToThreeD(models{i});
     317                                models{i}=TwoDToThreeD(models{i},self.planet);
    303318                        end
    304319                       
    … …  
    350365                function transfer(self,string) % {{{
    351366                        %Recover field size in one icecap:
    352                         eval(['n=length(self.icecaps{1}.' string ');']);
     367                        eval(['n=size(self.icecaps{1}.' string ',1);']);
    353368                        if n==self.icecaps{1}.mesh.numberofvertices,
    354369                                eval(['self.earth.' string '=zeros(self.earth.mesh.numberofvertices,1);']);
    … …  
    356371                                        eval(['self.earth.' string '(self.transitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
    357372                                end
     373                        elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
     374                                %dealing with a transient dataset.
     375                                %check that all timetags are similar between all icecaps:  %{{{
     376                                for i=1:length(self.icecaps),
     377                                        eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
     378                                        for j=(i+1):length(self.icecaps),
     379                                                eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
     380                                                if ~isequal(capfieldi(end,:),capfieldj(end,:)),
     381                                                        error(['Time stamps for ' string ' field is different between icecaps ' num2str(i) ' and ' num2str(j)]);
     382                                                end
     383                                        end
     384                                end
     385                                eval(['capfield1= self.icecaps{1}.' string ';']);
     386                                times=capfield1(end,:);
     387                                nsteps=length(times);
     388                                %}}}
     389                                %initialize:  %{{{
     390                                eval(['field=zeros(self.earth.mesh.numberofvertices+1,' num2str(nsteps) ');']);
     391                                eval(['field(end,:)=times;']); %transfer the times only, not the values
     392                                %}}}
     393                                %transfer all time fields: {{{
     394                                for i=1:length(self.icecaps),
     395                                        eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
     396                                        for j=1:nsteps,
     397                                                eval(['field(self.transitions{' num2str(i) '},' num2str(j) ')=capfieldi(1:end-1,' num2str(j) ');']); %transfer only the values, not the time.
     398                                        end
     399                                end
     400                                eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
     401                                %}}}
    358402                        elseif n==self.icecaps{1} .mesh.numberofelements,
    359403                                eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);

  • issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m

    r24155 r24779  
    1 function md=TwoDToThreeD(md)
     1function md=TwoDToThreeD(md,planet)
    22
    33        %reproject model into lat,long if necessary:
    4         if md.mesh.epsg~=4326,
    5                 [md.mesh.x,md.mesh.y]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),sprintf('EPSG:4326'));
     4        if ~strcmpi(md.mesh.proj,epsg2proj(4326)),
     5                [md.mesh.x,md.mesh.y]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
    66        end
    77
    88        %Make a 3dsurface mesh out of this:
    9         R=6.371012*10^6;
     9        R=planetradius(planet);
    1010
    1111        %we assume x and y hold the long,lat values:
    … …  
    1414        mask=md.mask;
    1515
     16        %Assume spherical body:
    1617        x = R .* cosd(lat) .* cosd(long);
    1718        y = R .* cosd(lat) .* sind(long);

  • issm/trunk-jpl/src/m/mesh/meshintersect3d.m

    r24766 r24779  
    6666        end
    6767
     68        if(~isempty(find(indices==0))) error('issue with transition vector having one empty slot'); end;
     69
    6870
    6971%               if length(s)>1,
Note: See TracChangeset for help on using the changeset viewer.