source: issm/trunk/src/m/classes/basin.m@ 24313

%BASIN class definition
%
%   Usage:
%      basin=basin();

classdef basin
        properties (SetAccess=public) 
                boundaries        = {};
                epsg              = 3426;
                name              = '';
                continent         = '';
        end
        methods (Static)
                function self = loadobj(self) % {{{
                        % This function is directly called by matlab when a model object is
                        % loaded. Update old properties here
                end% }}}
        end
        methods
                function self = basin(varargin) % {{{
                        switch nargin
                                case 0
                                        self=setdefaultparameters(self);
                                otherwise

                                        self=setdefaultparameters(self);
                                        options=pairoptions(varargin{:}); 
                        
                                        %recover field values: 
                                        self.boundaries=getfieldvalue(options,'boundaries',{});
                                        self.name=getfieldvalue(options,'name','');
                                        self.continent=getfieldvalue(options,'continent','');
                                        self.epsg=getfieldvalue(options,'epsg',3426);
                        end
                end % }}}
                function self = setdefaultparameters(self) % {{{
                        self.name='';
                        self.continent='';
                        self.epsg=3426;
                        self.boundaries={};

                end % }}}
                function disp(self) % {{{
                        disp(sprintf('   basin parameters:'));
                        fielddisplay(self,'continent','continent name');
                        fielddisplay(self,'name','basin name');
                        fielddisplay(self,'epsg','epsg projection number for the entire basin');
                        fielddisplay(self,'boundaries','list of boundary objects');
                        for i=1:length(self.boundaries),
                                disp(sprintf('             boundary #%i: %s',i,self.boundaries{i}.name));
                        end

                end % }}}
                function boolean=isnameany(self,varargin) % {{{
                        boolean=0;
                        for  i=1:length(varargin),
                                if strcmpi(self.name,varargin{i}), 
                                        boolean=1;
                                        break;
                                end
                        end
                end % }}}
                function boolean=iscontinentany(self,varargin) % {{{
                        boolean=0;
                        for  i=1:length(varargin),
                                if strcmpi(self.continent,varargin{i}), 
                                        boolean=1;
                                        break;
                                end
                        end
                end % }}}
                function output=outputname(self,varargin) % {{{
                
                        %recover options
                        options=pairoptions(varargin{:});
                        extension=getfieldvalue(options,'extension',1);

                        [path,nme,ext]=fileparts(self.name);
                        if extension,
                                output=[nme ext];
                        else
                                output=nme;
                        end
                end % }}}
                function plot(self,varargin) % {{{
        
                        %add option: 
                        for i=1:length(self.boundaries),
                                self.boundaries{i}.plot('epsg',self.epsg,varargin{:});
                        end

                end % }}}
                function plot3d(self,varargin) % {{{
        
                        %add option: 
                        for i=1:length(self.boundaries),
                                self.boundaries{i}.plot3d(varargin{:});
                        end

                end % }}}
                function out=contour(self,varargin) % {{{
                
                        %recover options
                        options=pairoptions(varargin{:});
                        x=[]; y=[];

                        %go through boundaries, recover edges and project them in the basin epsg reference frame: 
                        for i=1:length(self.boundaries),
                                boundary=self.boundaries{i};
                                contour=boundary.edges();
                                [contour.x,contour.y]=gdaltransform(contour.x,contour.y,sprintf('EPSG:%i',boundary.epsg),sprintf('EPSG:%i',self.epsg));
                                x=[x;contour.x];
                                y=[y;contour.y];
                        end
                        %close the contour: 
                        if x(end)~=x(1) | y(end)~=y(1), 
                                x(end)=x(1); y(end)=y(1);
                        end

                        out.x=x;
                        out.y=y;
                        out.nods=length(x);
                end % }}}
                function output=shapefilecrop(self,varargin) % {{{

                        %recover options
                        options=pairoptions(varargin{:});
                        threshold=getfieldvalue(options,'threshold',.65); %.65 degrees lat,long
                        inshapefile=getfieldvalue(options,'shapefile');
                        outputshapefile=getfieldvalue(options,'outputshapefile','');
                        epsgshapefile=getfieldvalue(options,'epsgshapefile');

                        %create list of contours that have critical length > threshold:  (in lat,long)
                        contours=shpread(inshapefile);
                        llist=[];
                        for i=1:length(contours),
                                contour=contours(i);
                                carea=polyarea(contour.x,contour.y);
                                clength=sqrt(carea);
                                if clength<threshold,
                                        llist=[llist;i];
                                end
                        end
                        contours(llist)=[];

                        %project onto reference frame:
                        if self.epsg~=epsgshapefile,
                                for i=1:length(contours),
                                        h=contours(i);
                                        [h.x,h.y]=gdaltransform(h.x,h.y,sprintf('EPSG:%i',epsgshapefile),sprintf('EPSG:%i',self.epsg));
                                        contours(i).x=h.x;
                                        contours(i).y=h.y;
                                end
                        end

                        %only keep the contours that are inside the basin (take centroids): 
                        ba=self.contour();
                        flags=zeros(length(contours),1);
                        for i=1:length(contours),
                                h=contours(i); 
                                in=inpolygon(h.x,h.y,ba.x,ba.y); 
                                if ~isempty(find(in==0)),
                                        flags(i)=1;
                                end
                        end
                        pos=find(flags);  contours(pos)=[];

                        %Two options: 
                        if strcmpi(outputshapefile,''),
                                output=contours;
                        else
                                shpwrite(contours,outputshapefile);
                        end

                end % }}}
        end
end