getfieldfromgrib_old.m

来自「麻省理工学院的人工智能工具箱,很珍贵,希望对大家有用!」· M 代码 · 共 140 行

function [TZH,dmn]=getFieldfromGRIB(ctl,dmn,varargin)

if(nargin<2 | isempty(dmn))
   [archivo,Camino2]=uigetfile('*.cfg','Selecciona un Dominio');
   dmn=readDomain([Camino2 archivo]);
end

startDate=datenum(dmn.startDate);
endDate=datenum(dmn.endDate);
if nargin>2,
   for i=1:2:length(varargin)
      switch lower(varargin{i}),
      case 'dates', 
         dates=varargin{i+1};
         if (datenum(dates{1})<startDate),
            warning(['no existen los datos correspondientes a la fecha ' dates{1}])
         else
            startDate = datenum(dates{1}); 
         end
         if (datenum(dates{2})>endDate),
            warning(['no existen los datos correspondientes a la fecha' dates{2}])
         else
            endDate= datenum(dates{2}); 
         end 
      end
   end
end


disp(dmn);
nP=size(dmn.par,2);
nH=size(dmn.tim,2);
nN=size(dmn.lvl,2);
%diamax=fix(dmn.tim(nH)/24);
diamax=max(fix(dmn.tim/24));

fcDate=datevec(startDate:endDate+diamax);
fcDate=fcDate(:,1:3);

ctlname=[ctl.cam ctl.fil]
fid=fopen(ctlname,'rb');
k=1;
while ~feof(fid)
   name{k}=[ctl.cam fgetl(fid)];
   k=k+1;
end
fclose(fid);
idxname=[ctlname '.idx'];
idxS=dir(idxname);
ctlS=dir(ctlname);


if (isempty(idxS) | datenum(idxS.date)<datenum(ctlS.date))
   buildIndex('make,sort,write',ctl.cam,ctl.fil,[ctl.fil '.idx']);   
else
   buildIndex('read',ctl.cam,ctl.fil,[ctl.fil '.idx']);
end

%buildIndex('read',ctl.cam,'','era.ctl.idx');

ndias=size(fcDate,1)-diamax;
%ndias=size(fcDate,1);
zona=dmn.nod;
nopun=size(zona,2);
disp(['Size(TZH)=' num2str([ndias nP*nH*nN*size(dmn.nod,2)])])
TZH(ndias,nP*nH*nN*size(dmn.nod,2))=0;
%DATA(size(fcDate,2),length(dmn.lon)*length(dmn.lat))=0;
size(TZH)
[xi,yi]=meshgrid(dmn.lon,dmn.lat);
for P=1:nP
   for H=1:nH
      hora=mod(dmn.tim(H),24);
      dia=fix(dmn.tim(H)/24);
      for N=1:nN
         disp([num2str([hora,dmn.par(P),0,dmn.lvl(N)])]);
         DATA=[];
         for D=1:size(fcDate,1)
            if(fcDate(D,2)==1 & fcDate(D,3)==1),fprintf(1,'%d\n',fcDate(D,1)),end
            mess=sprintf('%04d%02d%02d%02d%04d_%03d%03d%04d',fcDate(D,:),hora,0,dmn.par(P),0,dmn.lvl(N));
            F=buildIndex('find',mess,'1','');
            if F(1)
               [A,info]=readmessage(name{F(1)},F(2));
            else
               error(['No se ha encontrado un campo: ',mess]);
            end
            %if (info.PDS.Parameter==dmn.par(P) & info.PDS.Height1==dmn.lvl(N)& info.PDS.Hour==hora & info.PDS.Year==mod(fcDate(D,1),100) & info.PDS.Month==fcDate(D,2)& info.PDS.Day==fcDate(D,3))
            if (info.PDS.Parameter==dmn.par(P) & info.PDS.Height1==dmn.lvl(N)& info.PDS.Hour==hora & (info.PDS.Year+(info.PDS.Century-1)*100)==fcDate(D,1) & info.PDS.Month==fcDate(D,2)& info.PDS.Day==fcDate(D,3))
               %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               if isempty(DATA)
                  
                  scanmode=info.GDS.LatLon.ScanMode;
                  dx=(1-2*bitget(scanmode,1))*info.GDS.LatLon.Di;
                  dy=(-1+2*bitget(scanmode,2))*info.GDS.LatLon.Dj;
                  
                  x=info.GDS.LatLon.Lon1:dx:info.GDS.LatLon.Lon2;
                  y=info.GDS.LatLon.Lat1:dy:info.GDS.LatLon.Lat2;
                  x=x/1000;y=y/1000;
                  i=find(x>180);
                  x(i)=x(i)-360;
                  [X,Y]=meshgrid(x,y);   
                  if bitget(scanmode,3)
                     X=X';Y=Y';
                  end
                  DATA=zeros([size(fcDate,2),length(dmn.lon)*length(dmn.lat)])*NaN;
               end
               %  if (size(DATA,2) == prod(size(A)))
               %DATA(D,:)=A(:)';
               zi=interp2(X,Y,A,xi,yi,'linear')';
               DATA(D,:)=zi(:)';
               %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               %               A=A';
               %               if (size(DATA,2) == prod(size(A)))
               %                  DATA(D,:)=A(:)';
               %               else
               %                  error(sprintf('dim(A)=%d & dim(DATA,2)=%d are differents in mess: %s',prod(size(A)),size(DATA,2),mess));
               %               end
            else
              error('Problems reading message');
            end
            %DATA=load([dmn.src 'tim' num2str(hora,'%02d') '\par' num2str(dmn.par(P),'%03d') '_' num2str(dmn.lvl(N),'%04d')]);
         end
         
         %if (info.PDS.Parameter==dmn.par(kp) & info.PDS.Height1==dmn.lvl(kl) & info.PDS.PeriodTime==P1 & info.PDS.TimeInterval==P2)
         %   A=A';
         %   VC(kd,(1:nn)+nn*((kl-1)+nl*((kt-1)+nt*(kp-1))))=A(dmn.nod);
         %else
         %   error('Problems reading message');
         %end
         campo=DATA((1+dia):(end+dia-diamax),zona);
         %campo=DATA((1+dia):(end+dia-diamax),zona);
         %VC(1,(1:nn)+nn*(kl+nl*(kt+nt*kp)))=A(dmn.nod);
         ind=findVarPosition(dmn.par(P),dmn.tim(H),dmn.lvl(N),dmn);
         %ind=((P-1)*(nH*nN*nopun)+(H-1)*nN*nopun+(N-1)*nopun+1):((P-1)*(nH*nN*nopun)+(H-1)*nN*nopun+(N-1)*nopun+nopun);
         %ind=((H-1)*nN*nopun+(N-1)*nopun+1):((H-1)*nN*nopun+(N-1)*nopun+nopun);
         TZH(:,ind)=campo;
      end
   end
end