Matlab-处理OMI臭氧廓线数据

前段时间想用近地面的臭氧数据，站点观测数据年限太短，已有的卫星产品中大多是总柱浓度，找到了OMI的廓线数据，用的是近地层（一共18层，选择第1层，认为是近地层）

GES DISC 

 想记录的原因是该数据在处理时对我来讲费大劲了，首先1. 不像一般的nc数据在读取的时候直接通过ncdisp检查变量后输入变量名称即可，该nc数据是树状结构，在师兄的帮忙下用panoply看了一下数据结构

ncdisp下是这样

 

真正读取的时候按照它的结构这样即可读出，就是ncdisp中groups下的结构：

 2.这个数据在时间维度上是23，每天每53分钟一次记录，每次记录的是不同的经纬度下的臭氧数据（感觉有点类似于点数据），每幅文件里行数也不一致，和平常处理的格点数据不一样；采用插值（一天内所有记录）将该数据插为网格数据，再求月均值；

%%%
m=['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12'];
[xi,yi]=meshgrid(31:0.1:43,110:0.1:123);

for y=2005:2012

    for j=1:12
        filelist=dir(strcat('OMI-Aura_L2-OMO3PR_',num2str(y),'m',m(j,:),'*.nc4'));
        k=length(filelist);
        ozomm=zeros(131,121);
        index=zeros(131,121);
        for i=1:k
            filename=[filelist(i).name];
            ozo=ncread(filename,'/HDFEOS/SWATHS/O3Profile/Data Fields/O3');
            lat=ncread(filename,'/HDFEOS/SWATHS/O3Profile/Geolocation Fields/Latitude');
            lon=ncread(filename,'/HDFEOS/SWATHS/O3Profile/Geolocation Fields/Longitude');
            s=size(ozo);
            %%%S可能为空，if语句判断
            if length(s)<3
                index=index+zeros(131,121);
                ozomm=ozomm+zeros(131,121);
            else
                ozom=zeros(s(2),s(3));
                for ii=1:size(ozom,1)
                    for jj=1:size(ozom,2)
                        ozom(ii,jj)=ozo(:,ii,jj);
                    end
                end
                %%%插值处理
                ind=find(~isnan(ozom));
                [iii jjj]=ind2sub(size(ozom),ind);
                v=ozom(ind);
                latt=lat(ind);
                lonn=lon(ind);
                ozo2=griddata(latt,lonn,v,xi,yi);
                ozo2(find(isnan(ozo2)==1)) = 0;
                ozo2(find(ozo2<0)) = 0;
                %%%判断空集
                if isempty(ozo2)
                    index=index+zeros(131,121);
                    ozomm=ozomm+zeros(131,121);
                else
                    indd=ozo2~=0;
                    index=index+indd;
                    ozomm=ozomm+ozo2;
                end
            end
            
        end
        % 由于一些天数里部分像元缺失影像，应该除以有值的次数
        ozomm=ozomm./index;
        ozomm(find(isnan(ozomm)==1)) = 0;
        % 创建nc文件并导入经纬度信息
        InPath=strcat('OMI-Aura_L2-OMO3PR_',num2str(y),m(j,:),'.nc');
        nccreate(InPath,'lon','Dimensions',{'x',131});
        nccreate(InPath,'lat','Dimensions',{'y',121});
         % 读取经纬度信息
        lon2=yi(:,1);
        lat2=xi(1,:);
        ncwrite(InPath,'lon',lon2);
        ncwrite(InPath,'lat',lat2);
        % 创建变量并导入计算好的数据
        nccreate(InPath,'ozo',...
                 'Dimensions', {'x',131,'y',121},...
                 'FillValue','disable'); 
        ncwrite(InPath,'ozo',ozomm);
        %%%输出tif
        out='D:\ncc_datas\ideas\ERL\revise_data\o3\omi_o3\2004_2012\tif\';
        OutPath=strcat(out,'OMI-Aura_L3-OMSO2e_',num2str(y),m(j,:),'.tif');
        ozo3=ncread(InPath,'ozo');
        ozommm=rot90(ozo3,1);
        GeoRef = georasterref('Rastersize',size(ozommm),'Lonlim',[double(min(lon2)),double(max(lon2))],'Latlim',[double(min(lat2)),double(max(lat2))]);
        ozo_Tif =OutPath;
        geotiffwrite(ozo_Tif,flip(ozommm),GeoRef);
        
        %%%打印提示
        fprintf(strcat(num2str(y),m(j,:)))

    end


end

 处理完后和近地表臭氧数据集验证时，发现二者时间序列上不一致（坐标轴名称没有更换，实际为臭氧），在后续数值分析时结果不好；

但是还是记录一下这个费劲的过程~~~