经过这一段时间的对海洋数据的处理,接触了大量的与海洋相关的数据,例如海洋地形、海洋表面温度、盐度、湿度、云场、风场等数据,除了地形数据是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文件。
数据读取的代码如下:
头文件:
#ifndef CCMPFILEREAD_H
#define CCMPFILEREAD_H
class ccmpFileRead
{
public:
void ccmpFileRead::fileread(const char*ccmpFilename);
}; #endif // CCMPFILEREAD_H
源文件:
#include "ccmpfileread.h" #include <gdal_priv.h>
#include <vector>
#include <QVector> #include <string>
#include <QString>
#include <QStringList>
#include <QDebug> #include <fstream> using namespace std; void ccmpFileRead::fileread(const char *ccmpFilename)
{
vector <string> vFileSets;
vector <string> pStrDesc;
vector<vector<float>> allSSTPixelNum1,allSSTPixelNum2,allSSTPixelNum3; GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");//中文路径
GDALDataset* fileDataset = (GDALDataset*) GDALOpen(ccmpFilename,GA_ReadOnly);//打开HDF数据集
if (fileDataset == NULL)
{
return;
} char** sublist = GDALGetMetadata((GDALDatasetH) fileDataset,"SUBDATASETS");//获得数据的字符串,可以打印出来看看自己需要的数据在那 int iCount = CSLCount(sublist);
if(iCount <= ){
qDebug() << "该文件没有子数据" << endl;
GDALClose((GDALDriverH)fileDataset);
} //存储数据集信息
for(int i = ; sublist[i] != NULL;i++)
{ qDebug() << sublist[i] << endl; if(i% != )
{
continue;
} //三个数据集:uwnd vwnd wspd 只读取前两个数据集,第三个数据集是补充数据集 string tmpstr = sublist[i];
tmpstr = tmpstr.substr(tmpstr.find_first_of("=")+);
const char *tmpc_str = tmpstr.c_str(); string tmpdsc = sublist[i+];
tmpdsc = tmpdsc.substr(tmpdsc.find_first_of("=")+); GDALDataset* hTmpDt = (GDALDataset*)GDALOpen(tmpc_str,GA_ReadOnly);//打开该数据 if (hTmpDt != NULL)
{
vFileSets.push_back(tmpc_str);
}
if(&pStrDesc != NULL){
pStrDesc.push_back(tmpdsc);
}
GDALClose(hTmpDt);
} //三个数据集分别读取 qDebug() << "read uwnd ......" << endl; QString qtmpdsc1 = QString::fromStdString(pStrDesc[]);//锁定某一个数据集 qDebug()<<qtmpdsc1<<endl; float *lineData = NULL;
if (qtmpdsc1!=NULL)
{
GDALDataset *tempDt = (GDALDataset *)GDALOpen(vFileSets[].data(), GA_ReadOnly);
int BandNum = tempDt->GetRasterCount(); int panBandmap[] ={};
lineData = new float[ * *];
tempDt->RasterIO(GF_Read,,,,,lineData,,,GDT_Float32,,panBandmap,,,); for (int iLine = ; iLine <tempDt->GetRasterYSize(); iLine++)
{
allSSTPixelNum1.resize(tempDt->GetRasterYSize());
for (int iPixel = ; iPixel < tempDt->GetRasterXSize(); iPixel++)
{
allSSTPixelNum1[iLine].resize(tempDt->GetRasterXSize());
tempDt->RasterIO(GF_Read, , iLine, tempDt->GetRasterXSize(), ,lineData, tempDt->GetRasterXSize(), , GDT_Float32, , panBandmap,,,);
allSSTPixelNum1[iLine][iPixel] = lineData[iPixel];
} }
if(lineData)
{
delete[]lineData;
lineData = NULL;
} qDebug() << "uwnd read over!" << endl; qDebug() <<"uwnd="<<'\n'<<allSSTPixelNum1[]<<'\n'<<endl; } //d读取vwnd数据集 QString qtmpdsc2 = QString::fromStdString(pStrDesc[]); if (qtmpdsc2!=NULL)
{
GDALDataset *tempDt = (GDALDataset *)GDALOpen(vFileSets[].data(), GA_ReadOnly);
int BandNum = tempDt->GetRasterCount();
qDebug()<<BandNum<<endl;
int panBandmap[] ={};
lineData = new float[ * *];
tempDt->RasterIO(GF_Read,,,,,lineData,,,GDT_Float32,,panBandmap,,,); for (int iLine = ; iLine <tempDt->GetRasterYSize(); iLine++)
{
allSSTPixelNum2.resize(tempDt->GetRasterYSize());
for (int iPixel = ; iPixel < tempDt->GetRasterXSize(); iPixel++)
{
allSSTPixelNum2[iLine].resize(tempDt->GetRasterXSize());
tempDt->RasterIO(GF_Read, , iLine, tempDt->GetRasterXSize(), ,lineData, tempDt->GetRasterXSize(), , GDT_Float32, , panBandmap,,,);
allSSTPixelNum2[iLine][iPixel] = lineData[iPixel];
} }
if(lineData)
{
delete[]lineData;
lineData = NULL;
} qDebug() << "vwnd read over!" << endl; qDebug() <<"vwnd="<<'\n'<<allSSTPixelNum2[]<<'\n'<<endl; } //读取wspd数据 QString qtmpdsc3 = QString::fromStdString(pStrDesc[]); if (qtmpdsc3!=NULL)
{
GDALDataset *tempDt = (GDALDataset *)GDALOpen(vFileSets[].data(), GA_ReadOnly);
int BandNum = tempDt->GetRasterCount();
qDebug()<<BandNum<<endl;
int panBandmap[] ={};
lineData = new float[ * *];
tempDt->RasterIO(GF_Read,,,,,lineData,,,GDT_Float32,,panBandmap,,,); for (int iLine = ; iLine <tempDt->GetRasterYSize(); iLine++)
{
allSSTPixelNum3.resize(tempDt->GetRasterYSize());
for (int iPixel = ; iPixel < tempDt->GetRasterXSize(); iPixel++)
{
allSSTPixelNum3[iLine].resize(tempDt->GetRasterXSize());
tempDt->RasterIO(GF_Read, , iLine, tempDt->GetRasterXSize(), ,lineData, tempDt->GetRasterXSize(), , GDT_Float32, , panBandmap,,,);
allSSTPixelNum3[iLine][iPixel] = lineData[iPixel];
} } if(lineData)
{
delete[]lineData;
lineData = NULL;
} qDebug() << "wspd read over!" << endl; qDebug() <<"wspd="<<'\n'<<allSSTPixelNum3[]<<'\n'<<endl; GDALClose((GDALDatasetH)tempDt); } GDALClose((GDALDriverH)fileDataset);
}
主函数调用:
#include <QCoreApplication>
#include <ccmpfileread.h>
int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
ccmpFileRead a1;
a1.fileread("E:/odp_workplace/odp_data/testdata/CCMP_Wind_Analysis_198707_V02.0_L3.5_RSS.nc");
return a.exec();
}
输出结果:
如上图所示数据已经读取并显示成功。