ps_plot.m

来自「StaMps最新测试版」· M 代码 · 共 796 行 · 第 1/2 页

        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; % m(1,:) is master APS + mean deviation from model
        %ph_all=m(3,:)'/4/pi*lambda*1000;
        textsize=0;
        units='mm/yr';
     case {'v-da'}
        uw=load(phuwname);
        scla=load(sclaname);
        aps=load(apsname);
        ph_uw=uw.ph_uw - scla.ph_scla- aps.ph_aps_slave;
        clear uw scla aps
        ph_all=zeros(n_ps,1);
        ref_ps=ps_setref;
        unwrap_ifg_index=setdiff(unwrap_ifg_index,ps.master_ix);
        if ~isempty(ifg_list)
            unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
            ifg_list=[];
        end
        ph_uw=ph_uw(:,unwrap_ifg_index);
        day=day(unwrap_ifg_index);
        
        ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),n_ps,1);
        %G=[ones(size(day)),day-master_day,cos(2*pi*(day-datenum(2000,1,0))/365.25-5.86)]; 
        G=[ones(size(day)),day-master_day]; 
        
        m=G\double(ph_uw');
        save mean_v m
        lambda=getparm('lambda');
        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; % m(1,:) is master APS + mean deviation from model
        %ph_all=m(3,:)'/4/pi*lambda*1000;
        textsize=0;
        units='mm/yr';
     case {'v-a'}
        uw=load(phuwname);
        aps=load(apsname);
        ph_uw=uw.ph_uw - aps.ph_aps_slave;
        clear uw aps
        ph_all=zeros(n_ps,1);
        ref_ps=ps_setref;
        unwrap_ifg_index=setdiff(unwrap_ifg_index,ps.master_ix);
        if ~isempty(ifg_list)
            unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
            ifg_list=[];
        end
        ph_uw=ph_uw(:,unwrap_ifg_index);
        day=day(unwrap_ifg_index);
        
        ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),n_ps,1);
        %G=[ones(size(day)),day-master_day,cos(2*pi*(day-datenum(2000,1,0))/365.25-5.86)]; 
        G=[ones(size(day)),day-master_day]; 
        
        m=G\double(ph_uw');
        save mean_v m
        lambda=getparm('lambda');
        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; % m(1,:) is master APS + mean deviation from model
        %ph_all=m(3,:)'/4/pi*lambda*1000;
        textsize=0;
        units='mm/yr';
     case {'v-ds'}
        uw=load(phuwname);
        scla=load(sclaname);
        scn=load(scnname);
        ph_uw=uw.ph_uw - scla.ph_scla - scn.ph_scn_slave;
        clear uw scla scn

        ph_all=zeros(n_ps,1);
        ref_ps=ps_setref;
        unwrap_ifg_index=setdiff(unwrap_ifg_index,ps.master_ix);
        if ~isempty(ifg_list)
            unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
            ifg_list=[];
        end
        ph_uw=ph_uw(:,unwrap_ifg_index);
        day=day(unwrap_ifg_index);
        
        ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),n_ps,1);
        %G=[ones(size(day)),day-master_day,cos(2*pi*(day-datenum(2000,1,0))/365.25-5.86)]; 
        G=[ones(size(day)),day-master_day]; 
        
        m=G\double(ph_uw');
        save mean_v m
        lambda=getparm('lambda');
        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; % m(1,:) is master APS + mean deviation from model
        %ph_all=m(3,:)'/4/pi*lambda*1000;
        textsize=0;
        units='mm/yr';
     case {'v-das'}
        uw=load(phuwname);
        scla=load(sclaname);
        scn=load(scnname);
        aps=load(apsname);
        ph_uw=uw.ph_uw - scla.ph_scla- aps.ph_aps_slave - scn.ph_scn_slave;
        clear uw scla aps scn

        ph_all=zeros(n_ps,1);
        ref_ps=ps_setref;
        unwrap_ifg_index=setdiff(unwrap_ifg_index,ps.master_ix);
        if ~isempty(ifg_list)
            unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
            ifg_list=[];
        end
        ph_uw=ph_uw(:,unwrap_ifg_index);
        day=day(unwrap_ifg_index);
        
        ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),n_ps,1);
        %G=[ones(size(day)),day-master_day,cos(2*pi*(day-datenum(2000,1,0))/365.25-5.86)]; 
        G=[ones(size(day)),day-master_day]; 
        
        m=G\double(ph_uw');
        save mean_v m
        lambda=getparm('lambda');
        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; % m(1,:) is master APS + mean deviation from model
        %ph_all=m(3,:)'/4/pi*lambda*1000;
        textsize=0;
        units='mm/yr';
   case {'vdrop-d'}
        uw=load(phuwname);
        scla=load(sclaname);
        ph_uw=uw.ph_uw - scla.ph_scla;
        clear uw scla
        unwrap_ifg_index=setdiff(unwrap_ifg_index,ps.master_ix);
        ph_all=zeros(n_ps,1);
        ref_ps=ps_setref;
        if ~isempty(ifg_list)
            unwrap_ifg_index=intersect(unwrap_ifg_index,ifg_list);
            ifg_list=[];
        end
        ph_uw=ph_uw(:,unwrap_ifg_index);
        day=day(unwrap_ifg_index);
        
        ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),n_ps,1);
        %G=[ones(size(day)),day-master_day,cos(2*pi*(day-datenum(2000,1,0))/365.25-5.86)]; 
        G=[ones(size(day)),day-master_day] ;
        n=size(ph_uw,2);
        lambda=getparm('lambda');
        ph_all=zeros(size(ph_uw));
        for i=1:n
            m=G([1:i-1,i+1:end],:)\double(ph_uw(:,[1:i-1,i+1:n])');
            ph_all(:,i)=m(2,:)'*365.25/4/pi*lambda*1000; 
        end
        %ph_all=m(3,:)'/4/pi*lambda*1000;
        textsize=0;     
        units='mm/yr';
   case {'vsb'}
        phuw=load(phuwsbname);
        ref_ps=ps_setref;
        unwrap_ifg_index=getparm('unwrap_ifg_index')
        if strcmp(unwrap_ifg_index,'all')
            unwrap_ifg_index=[1:n_ifg];
        end	
        ph_uw=phuw.ph_uw(:,unwrap_ifg_index);
        ifgday_ix=ps.ifgday_ix(unwrap_ifg_index,:);
        clear phuw
        ph_uw=ph_uw-repmat(mean(ph_uw(ref_ps,:)),n_ps,1);
        G=[ones(size(ifgday_ix(:,1))),day(ifgday_ix(:,2))-day(ifgday_ix(:,1))];
        m=G\double(ph_uw');
        lambda=getparm('lambda');
        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; 
        %ph_all=m(1,:)'; 
        textsize=0;
        units='mm/yr';
   case {'vsb-d'}
        uw=load(phuwsbname);
        scla=load(sclasbname);
        ph_all=uw.ph_uw - scla.ph_scla;
        clear uw scla
        ref_ps=ps_setref;
        unwrap_ifg_index=getparm('unwrap_ifg_index')
        if strcmp(unwrap_ifg_index,'all')
            unwrap_ifg_index=[1:n_ifg];
        end	
        ph_all=ph_all(:,unwrap_ifg_index);
        ifgday_ix=ps.ifgday_ix(unwrap_ifg_index,:);
        clear phuw
        ph_all=ph_all-repmat(mean(ph_all(ref_ps,:)),n_ps,1);
        G=[ones(size(ifgday_ix(:,1))),day(ifgday_ix(:,2))-day(ifgday_ix(:,1))];
        m=G\double(ph_all');
        lambda=getparm('lambda');
        ph_all=m(2,:)'*365.25/4/pi*lambda*1000; 
        %ph_all=m(1,:)'; 
        textsize=0;
        units='mm/yr';
    case {'p'}
        pm=load(pmname);
        %pm=load(rcuwname);
        ph_all=pm.ph_patch./abs(pm.ph_patch);
        if n_ifg~=size(ph_all,2)
            ph_all=[ph_all(:,1:ps.master_ix-1),zeros(ps.n_ps,1),ph_all(:,ps.master_ix:end)];
        end
        clear pm
        if ref_ifg~=0
            ph_all=ph_all.*repmat(conj(ph_all(:,ref_ifg)),1,n_ifg);
            ph_all(:,ref_ifg)=1; % may not be so because of rounding errors
        end
    case {'wf'}
        uw=load('uw_grid');
        gridix=zeros(size(uw.nzix));
        gridix(uw.nzix)=[1:uw.n_ps];
        ph_all=zeros(ps.n_ps,uw.n_ifg);
        for i=1:ps.n_ps
            ph_all(i,:)=uw.ph(gridix(uw.grid_ij(i,1),uw.grid_ij(i,2)),:);
        end
        clear uw
        if ref_ifg~=0
            ph_all=ph_all.*conj(repmat(ph_all(:,ref_ifg),1,n_ifg));
        end
    case {'vs-d'}
        if exist(['./',meanvname,'.mat'],'file')==0
            ps_mean_v;
        end
        mv=load(meanvname);
        ph_all=mv.mean_v_std;
        units='mm/yr';

        
        
    otherwise
        error('unknown value type')
end

if isempty(ifg_list)
    ifg_list=1:size(ph_all,2);
end
n_ifg_plot=length(ifg_list);

figure
set(gcf,'renderer','zbuffer')

[Y,X]=meshgrid([0.7:-0.2:0.1],[0.1:0.1:0.8]);
xgap=0.1;
ygap=0.2;

if ~isempty(lon_rg)
    ix=lonlat(:,1)>=lon_rg(1)&lonlat(:,1)<=lon_rg(2);
    lonlat=lonlat(ix,:);
end

if ~isempty(lat_rg)
    ix=lonlat(:,2)>=lat_rg(1)&lonlat(:,2)<=lat_rg(2);
    lonlat=lonlat(ix,:);
end

max_xy=llh2local([max(lonlat),0]',[min(lonlat),0]);
fig_ar=4/3; % aspect ratio of figure window
useratio=1; % max fraction of figure window to use
n_i=max_xy(2)*1000;
n_j=max_xy(1)*1000;
ar=max_xy(1)/max_xy(2); % aspect ratio (x/y)
if n_x==0
    n_y=ceil(sqrt((n_ifg_plot)*ar/fig_ar)); % number of plots in y direction
    n_x=ceil((n_ifg_plot)/n_y);
else
    n_y=ceil((n_ifg_plot)/n_x);
end


% d_y=0.9/n_y;
% d_x=d_y*ar;
% if d_x>0.9/n_x
%     d_x=0.9/n_x;
%     d_y=d_x/ar;
% end
d_x=useratio/n_x;
d_y=d_x/ar*fig_ar;
if d_y>useratio/n_y
    d_y=useratio/n_y; 
    d_x=d_y*ar/fig_ar;
end

h_y=0.95*d_y;
h_x=h_y*ar/fig_ar;
%y=1-d_y:-0.99/n_y:0;
%y=1-d_y:-0.90/n_y:0;
y=1-d_y:-d_y:0;
x=1-useratio:d_x:1-d_x;


[imY,imX]=meshgrid(y,x);
if textsize==0
    textsize=round(10*4/n_x);
    if textsize>16
        textsize=16;
    elseif textsize<8
        textsize=8;
    end
end

l_t=1/9*abs(textsize)/10; % text length
h_t=1/50*abs(textsize)/10; % text height
x_t=round((h_x-l_t)/h_x/2*n_j);
y_t=round(h_t*1.2/h_y*n_i);



ph_disp=ph_all(:,ifg_list);

if isreal(ph_all)

    if ref_ifg~=0
        if ref_ifg==-1
            ph_disp=ph_disp-[ph_disp(:,1),ph_disp(:,1:end-1)];
        else
            ph_disp=ph_disp-repmat(ph_all(:,ref_ifg),1,size(ph_disp,2));
        end
    else
        ref_ifg=master_ix;
    end
    if ref_ps~=0
        ref_ph=(ph_disp(ref_ps,:));
        mean_ph=zeros(1,size(ph_disp,2));
        for i=1:size(ph_disp,2)
            mean_ph(i)=mean(ref_ph(~isnan(ref_ph(:,i)),i));
        end
        ph_disp=ph_disp-repmat(mean_ph,n_ps,1);
    end
    
    phsort=sort(ph_disp(~isnan(ph_disp)));
    if isempty(lims)
        maxph=phsort(round(length(phsort)*.999));
        minph=phsort(ceil(length(phsort)*.001));
       lims=[minph,maxph];
    end
else
    if ref_ifg==0
        ref_ifg=master_ix;
    elseif ref_ifg==-1
        ph_disp=ph_disp.*conj([ph_disp(:,1),ph_disp(:,1:end-1)]);
    end
    if ref_ps~=0
        ph_disp=ph_disp./abs(ph_disp);
        ref_ph=(ph_disp(ref_ps,:));
        mean_ph=zeros(1,size(ph_disp,2));
        for i=1:size(ph_disp,2)
            mean_ph(i)=sum(ref_ph(~isnan(ref_ph(:,i)),i));
        end
        ph_disp=ph_disp.*conj(repmat(mean_ph,n_ps,1));
    end
    lims=[-pi,pi];

end

%for i=1:size(ph_all,2)
i_im=0;
for i=ifg_list
    %subplot(5,7,i);
    i_im=i_im+1;
    if n_ifg_plot>1
        axes('position',[imX(i_im),imY(i_im),h_x,h_y])
    end

    %axes('position',[X(i),Y(i),xgap,ygap])
    ps_plot_ifg(ph_disp(:,i_im),plot_flag,lims,lon_rg,lat_rg);
    %plot_phase(ph_tc(:,i)*conj(ph_tc(ref_ix,i)));
    box on
    if n_ifg_plot>1
        set(gca,'yticklabel',[])
        set(gca,'xticklabel',[])
    end
    xlim=get(gca,'xlim');
    x_t=(h_x-l_t)/2/h_x*(xlim(2)-xlim(1))+xlim(1);
    ylim=get(gca,'ylim');
    if textsize>0
        y_t=(h_y-1.2*h_t)/h_y*(ylim(2)-ylim(1))+ylim(1);
    else
        y_t=(0.5*h_t)/h_y*(ylim(2)-ylim(1))+ylim(1);
    end
    %xlabel([num2str((day(i)/365.25),3),'yr, ',num2str(round(bperp(i))),'m'])
    if textsize~=0 & size(day,1)==size(ph_all,2)
        t=text(x_t,y_t,[datestr(day(i),'dd mmm yyyy')]);
        set(t,'fontweight','bold','color',textcolor,'fontsize',abs(textsize))
    end
    
    if cbar_flag==0 & (i==ref_ifg | (isempty(intersect(ref_ifg,ifg_list)) & i==ifg_list(1))) 
        if n_ifg_plot>1
            h=colorbar('South');
        else
            h=colorbar('SouthOutside');
        end
        xlim=get(h,'xlim');
        set(h,'xlim',[xlim(2)-64,xlim(2)])
        if diff(lims)>1 | diff(lims)==0
            plotlims=round(lims*10)/10;
        else
            limorder=ceil(-log10(diff(lims)))+2;
            plotlims=round(lims*10^limorder)/10^limorder;
        end
        set(h,'xtick',[xlim(2)-64,xlim(2)],'Xticklabel',plotlims,'xcolor',textcolor,'ycolor',textcolor,'fontweight','bold','color',textcolor)
    end

end
    

fprintf('Color Range: %g to %g %s\n',lims,units)