欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

基于GDAL库,读取海洋风场数据(.nc格式)c++版

程序员文章站 2022-05-30 21:59:00
经过这一段时间的对海洋数据的处理,接触了大量的与海洋相关的数据,例如海洋地形、海洋表面温度、盐度、湿度、云场、风场等数据,除了地形数据是grd格式外,其他的都是nc格式的数据。本文将以海洋风场数据为例,进行nc格式文件的读取。 海洋风场数据(ccmp_wind)一般情况下会包含三个数据集:第一个数据 ......

        经过这一段时间的对海洋数据的处理,接触了大量的与海洋相关的数据,例如海洋地形、海洋表面温度、盐度、湿度、云场、风场等数据,除了地形数据是grd格式外,其他的都是nc格式的数据。本文将以海洋风场数据为例,进行nc格式文件的读取。

         海洋风场数据(ccmp_wind)一般情况下会包含三个数据集:第一个数据集是uwnd(standard_name = "eastward_wind"),第二个数据集是vwnd(standard_name = "northward_wind"),第三个数据集是nobs或者wspd。前两个数据集是矢量数据,表示此处的风场方向最后一个数据集是标量数据,代表此处的风速。每个数据集中数据的存储又分为四个波段(也可以说是图层),一天的观测时间分为四个时间点,所以有四个图层。

          gdal库可以提供对nc格式数据的读取,本次数据的读取是在qt+vs2017环境下配置gdal库和netcdf库,环境的配置可以在网上找到,gdal库的配置可以根据《gdal源码剖析和开发指南》书中的内容进行编译和配置,配置完成后就可以运行数据,读取nc文件。

           数据读取的代码如下:

头文件:

 1 #ifndef ccmpfileread_h
 2 #define ccmpfileread_h
 3 class ccmpfileread
 4 {
 5 public:
 6     void ccmpfileread::fileread(const char*ccmpfilename);
 7 };
 8 
 9 
10 
11 #endif // ccmpfileread_h

源文件:

  1 #include "ccmpfileread.h"
  2 
  3 #include <gdal_priv.h>
  4 #include <vector>
  5 #include <qvector>
  6 
  7 #include <string>
  8 #include <qstring>
  9 #include <qstringlist>
 10 #include <qdebug>
 11 
 12 #include <fstream>
 13 
 14 using namespace std;
 15 
 16 void ccmpfileread::fileread(const char *ccmpfilename)
 17 {
 18     vector <string>         vfilesets;
 19     vector <string>         pstrdesc;
 20     vector<vector<float>>   allsstpixelnum1,allsstpixelnum2,allsstpixelnum3;
 21 
 22 
 23     gdalallregister();
 24     cplsetconfigoption("gdal_filename_is_utf8","no");//中文路径
 25     gdaldataset* filedataset = (gdaldataset*) gdalopen(ccmpfilename,ga_readonly);//打开hdf数据集
 26     if (filedataset == null)
 27     {
 28         return;
 29     }
 30 
 31     char** sublist = gdalgetmetadata((gdaldataseth) filedataset,"subdatasets");//获得数据的字符串,可以打印出来看看自己需要的数据在那
 32 
 33     int icount = cslcount(sublist);
 34     if(icount <= 0){
 35         qdebug() << "该文件没有子数据" << endl;
 36         gdalclose((gdaldriverh)filedataset);
 37     }
 38 
 39     //存储数据集信息
 40     for(int i = 0; sublist[i] != null;i++)
 41     {
 42 
 43         qdebug() << sublist[i] << endl;
 44 
 45         if(i%2 != 0)
 46         {
 47             continue;
 48         }
 49 
 50         //三个数据集:uwnd vwnd wspd 只读取前两个数据集,第三个数据集是补充数据集
 51 
 52         string tmpstr = sublist[i];
 53         tmpstr = tmpstr.substr(tmpstr.find_first_of("=")+1);
 54         const char *tmpc_str = tmpstr.c_str();
 55 
 56         string tmpdsc = sublist[i+1];
 57         tmpdsc = tmpdsc.substr(tmpdsc.find_first_of("=")+1);
 58 
 59         gdaldataset* htmpdt = (gdaldataset*)gdalopen(tmpc_str,ga_readonly);//打开该数据
 60 
 61         if (htmpdt != null)
 62         {
 63             vfilesets.push_back(tmpc_str);
 64         }
 65         if(&pstrdesc != null){
 66             pstrdesc.push_back(tmpdsc);
 67         }
 68         gdalclose(htmpdt);
 69     }
 70 
 71 
 72 //三个数据集分别读取
 73 
 74     qdebug() << "read uwnd ......" << endl;
 75 
 76     qstring qtmpdsc1 = qstring::fromstdstring(pstrdesc[0]);//锁定某一个数据集
 77 
 78     qdebug()<<qtmpdsc1<<endl;
 79 
 80     float *linedata = null;
 81     if (qtmpdsc1!=null)
 82     {
 83         gdaldataset  *tempdt = (gdaldataset *)gdalopen(vfilesets[0].data(), ga_readonly);
 84                int bandnum = tempdt->getrastercount();
 85 
 86                int panbandmap[1] ={1};
 87                linedata = new float[1 * 200*200];
 88           tempdt->rasterio(gf_read,508,112,32,24,linedata,50,50,gdt_float32,1,panbandmap,0,0,0);
 89 
 90 
 91            for (int iline = 0; iline <tempdt->getrasterysize(); iline++)
 92             {
 93                  allsstpixelnum1.resize(tempdt->getrasterysize());
 94             for (int ipixel = 0; ipixel < tempdt->getrasterxsize(); ipixel++)
 95              {
 96                   allsstpixelnum1[iline].resize(tempdt->getrasterxsize());
 97                   tempdt->rasterio(gf_read, 0, iline, tempdt->getrasterxsize(), 1,linedata, tempdt->getrasterxsize(), 1, gdt_float32, 1, panbandmap,0,0,0);
 98                   allsstpixelnum1[iline][ipixel] = linedata[ipixel];
 99                 }
100 
101            }
102            if(linedata)
103              {
104             delete[]linedata;
105             linedata = null;
106               }
107 
108            qdebug() << "uwnd read over!" << endl;
109 
110            qdebug() <<"uwnd="<<'\n'<<allsstpixelnum1[200]<<'\n'<<endl;
111 
112         }
113 
114     //d读取vwnd数据集
115 
116     qstring qtmpdsc2 = qstring::fromstdstring(pstrdesc[2]);
117 
118     if (qtmpdsc2!=null)
119     {
120         gdaldataset  *tempdt = (gdaldataset *)gdalopen(vfilesets[0].data(), ga_readonly);
121                int bandnum = tempdt->getrastercount();
122                qdebug()<<bandnum<<endl;
123                int panbandmap[1] ={1};
124                linedata = new float[1 * 200*200];
125           tempdt->rasterio(gf_read,508,112,32,24,linedata,50,50,gdt_float32,1,panbandmap,0,0,0);
126 
127 
128            for (int iline = 0; iline <tempdt->getrasterysize(); iline++)
129             {
130                  allsstpixelnum2.resize(tempdt->getrasterysize());
131             for (int ipixel = 0; ipixel < tempdt->getrasterxsize(); ipixel++)
132              {
133                   allsstpixelnum2[iline].resize(tempdt->getrasterxsize());
134                   tempdt->rasterio(gf_read, 0, iline, tempdt->getrasterxsize(), 1,linedata, tempdt->getrasterxsize(), 1, gdt_float32, 1, panbandmap,0,0,0);
135                   allsstpixelnum2[iline][ipixel] = linedata[ipixel];
136                 }
137 
138            }
139            if(linedata)
140              {
141             delete[]linedata;
142             linedata = null;
143               }
144 
145            qdebug() << "vwnd read over!" << endl;
146 
147            qdebug() <<"vwnd="<<'\n'<<allsstpixelnum2[200]<<'\n'<<endl;
148 
149         }
150 
151     //读取wspd数据
152 
153     qstring qtmpdsc3 = qstring::fromstdstring(pstrdesc[2]);
154 
155     if (qtmpdsc3!=null)
156     {
157         gdaldataset  *tempdt = (gdaldataset *)gdalopen(vfilesets[0].data(), ga_readonly);
158                int bandnum = tempdt->getrastercount();
159                qdebug()<<bandnum<<endl;
160                int panbandmap[1] ={1};
161                linedata = new float[1 * 200*200];
162           tempdt->rasterio(gf_read,508,112,32,24,linedata,50,50,gdt_float32,1,panbandmap,0,0,0);
163 
164 
165            for (int iline = 0; iline <tempdt->getrasterysize(); iline++)
166             {
167                  allsstpixelnum3.resize(tempdt->getrasterysize());
168             for (int ipixel = 0; ipixel < tempdt->getrasterxsize(); ipixel++)
169              {
170                   allsstpixelnum3[iline].resize(tempdt->getrasterxsize());
171                   tempdt->rasterio(gf_read, 0, iline, tempdt->getrasterxsize(), 1,linedata, tempdt->getrasterxsize(), 1, gdt_float32, 1, panbandmap,0,0,0);
172                   allsstpixelnum3[iline][ipixel] = linedata[ipixel];
173                 }
174 
175            }
176 
177            if(linedata)
178              {
179             delete[]linedata;
180             linedata = null;
181               }
182 
183            qdebug() << "wspd read over!" << endl;
184 
185            qdebug() <<"wspd="<<'\n'<<allsstpixelnum3[200]<<'\n'<<endl;
186 
187            gdalclose((gdaldataseth)tempdt);
188 
189         }
190 
191         gdalclose((gdaldriverh)filedataset);
192 }

主函数调用:

1 #include <qcoreapplication>
2 #include <ccmpfileread.h>
3 int main(int argc, char *argv[])
4 {
5     qcoreapplication a(argc, argv);
6     ccmpfileread a1;
7     a1.fileread("e:/odp_workplace/odp_data/testdata/ccmp_wind_analysis_198707_v02.0_l3.5_rss.nc");
8     return a.exec();
9 }

输出结果:

基于GDAL库,读取海洋风场数据(.nc格式)c++版基于GDAL库,读取海洋风场数据(.nc格式)c++版基于GDAL库,读取海洋风场数据(.nc格式)c++版

如上图所示数据已经读取并显示成功。