利用Matlab编程实现对NetCDF数据针对某一坐标按时间提取风、流、水蒸气等变量信息

一、NetCDF数据介绍

       NetCDF(network Common Data Form)网络通用数据格式是由美国大学大气研究协会（University Corporation for Atmospheric Research，UCAR)的Unidata项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。用户可以借助多种方式方便地管理和操作 NetCDF 数据集。

二、数据下载

1. 进入网址 https://www.ecmwf.int/ 申请一个账号，注册流程简单，按操作提示完成即可。
2. 文件下载，找到数据集，填写检索参数，即可下载。               

                                          

三、数据提取需求

利用Matlab编程实现提取**地区三个坐标点的风（u10、v10）、总降雨(tp)、水蒸气含量(tcwv)信息，且输出结果组织如下表：

提取结果组织形式

时间	u10	v10	tp	tcwv
20100101
20100102
……
20130101
……
……

四、数据提取思想

要想利用编程提取，我们首先需要学习一下NetCDF的文件组织形式：

1、读取NetCDF的源数据信息。在Matlab2015带有读取该文件类型的函数:ncdisp('path')，示例如下：

ncdisp('C:\Users\Administrator\Desktop\dataset\2017\1月份\_grib2netcdf-webmars-public-svc-green-002-6fe5cac1a363ec1525f54343b6cc9fd8-4vMyn2.nc') 

2、读取NetCDF的变量信息。在Matlab2015带有读取该文件变量的函数：ncread('path','变量的名称')，示例如下： 

ncread('C:\Users\Administrator\Desktop\dataset\2017\1月份\_grib2netcdf-webmars-public-svc-green-002-6fe5cac1a363ec1525f54343b6cc9fd8-4vMyn2.nc','time')
ncread('C:\Users\Administrator\Desktop\dataset\2017\1月份\_grib2netcdf-webmars-public-svc-green-002-6fe5cac1a363ec1525f54343b6cc9fd8-4vMyn2.nc','longitude')

3、由上面提取的变量，在Matlab查看可知：

          变量数据结构：由下知，u10、v10、tp、tcwv变量空间索引为    经度   维度   时间
          time           30x1 int32
          latitude      721x1 single （可知经纬度间隔，方便推算对应经纬度下变量的空间坐标索引）
          longitude   1440x1 single（可知经纬度间隔，方便推算对应经纬度下变量的空间坐标索引）
          u10            1440x721x30 double
          v10            1440x721x30 double
          tp              1440x721x30 double
          tcwv           1440x721x30 double

五、Matlab代码实现

clear 
close all

savedir = 'I:\feng_liu_data\result\';  % 保存excel的路径
point = [50 3;93 1.5;120 6.5];     % 经纬度的坐标点
startYear = 2017;endYear = 2019;    % 数据起始年份

%遍历文件夹找寻下载的NetCDF文件
j = 1;
for year = startYear : endYear  % 遍历年份
    maindir = ['I:\feng_liu_data\datset\' num2str(year) '\'];   % nc文件的路径，是年份再上一级的路径，因为我的路径是：I:\feng_liu_data\datset\里面有三个文件夹（2017、2018、2019），每个年份文件夹里有12个月文件夹，每个文件夹里才是NC文件
    subdir =  dir( maindir );      % 读取maindir文件夹中的所有文件名，返回一列
    for i = 1 : length( subdir )
        if( isequal( subdir( i ).name, '.' ) || ...
            isequal( subdir( i ).name, '..' ) || ...
            ~subdir( i ).isdir )   % 如果不是目录跳过
            continue;
        end
        subdirpath1 = fullfile( maindir, subdir( i ).name, '*.nc' );
        file = dir( subdirpath1 );   % 在这个子文件夹下找后缀为nc的文件    
        for jj = 1 : length( file )  % 遍历文件
            filedir{j} = fullfile( maindir, subdir( i ).name);  % 保存文件目录
            filepath{j} = fullfile( maindir, subdir( i ).name, file( jj ).name  );  % 保存文件路径
            filename{j} =  file(jj).name;       % 保存文件名字
            j = j + 1;
        end
    end 
end


% 建立一些空的变量，用来记录所有时间下的数据
year_u10 = [];year_v10 = [];year_tp = [];year_tcwv = [];year_time = [];month = 1;thisYear = startYear;
for i = 1 : j-1
    % 读取NC文件的时间、经纬度
    time = ncread(filepath{i},'time'); 
    longitude = ncread(filepath{i},'longitude');
    latitude = ncread(filepath{i},'latitude');
    % 找到对应坐标下的序号
    for k = 1 : size(point,1)
        index_lon(k) = (point(k,1)/0.25)+1;
        index_lat(k) = (90-point(k,2))/0.25+1;
    end
    % 读取要保存的信息
    u10 = ncread(filepath{i},'u10');v10 = ncread(filepath{i},'v10');tp = ncread(filepath{i},'tp');tcwv = ncread(filepath{i},'tcwv');
    % 将对应经纬度下的所有时间变量数据筛选出存入下面变量中
    new_u10 =[];new_v10 =[];new_tp =[];new_tcwv =[];
    for k = 1 : size(point,1)
        new_u10(k,:) = u10(index_lon(k),index_lat(k),:);
        new_v10(k,:) = v10(index_lon(k),index_lat(k),:);
        new_tp(k,:) = tp(index_lon(k),index_lat(k),:);
        new_tcwv(k,:) = tcwv(index_lon(k),index_lat(k),:);
    end
    % 写入到刚刚这些变量，用于最后输出
    year_u10 = [year_u10 new_u10];year_v10 = [year_v10 new_v10];year_tp = [year_tp new_tp];year_tcwv = [year_tcwv new_tcwv];
    % 记录时间
    n = length(time);
     thisTime = num2str(thisYear*10000 + month*ones(n,1)*100 + [1:n]');%根据已有数据是2017年到2019年，按照time个数推算日期
     year_time = [year_time;thisTime];
    %计算月份
    month = month + 1;
    if month == 13   %判断一年是否结束
        thisYear = thisYear + 1;
        month = 1;
    end
end
   q=[];r=[];%用来存放输出EXCLE的名称的坐标点
for k = 1 : size(point,1)
    data = cell( length(year_time)+1,5);%由于存在文本和数字，建立cell元胞数组
    data{1,1} = '时间';data{1,2} = 'u10';data{1,3} = 'V10';data{1,4} = 'tp';data{1,5} = 'tcwv';%表头存入元胞数组第一行中
    for i = 1 : length(year_time)
        data(i+1,:) = {year_time(i,:) year_u10(k,i)' year_v10(k,i)' year_tp(k,i)' year_tcwv(k,i)'};
        q = point(k);
        r = point(3+k);
    end
    
    xlswrite(['I:\feng_liu_data\result\' 'point' '(' num2str(q) ',' num2str(r) ')' '.xlsx'],data)
end

六、提取结果

                              

七、总结

静下心认真思考，相信你可以。