基于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 }
输出结果:
如上图所示数据已经读取并显示成功。
上一篇: 拉格朗日插值学习小结