C++读取NC数据的结果与真实数值不一致的解决方法

  本文介绍基于C++ 语言的netCDF库读取.nc格式的栅格文件时，代码读取到的数据与栅格文件的实际数据不一致的解决方法。

  最近，由于需要读取ERA5气象数据，因此使用C++ 语言中的netCDF库读取.nc格式文件。其中，偶然发现在Visual Studio的代码中读取到的.nc文件的数据，和其实际的数据（ArcMap等软件打开所显示的数据）不一致；这里就介绍一种可能导致上述情况的原因，以及对应的解决方法。此外，关于Visual Studio中配置C++ 语言netCDF库的方法，大家可以参考在Visual Studio配置C++的netCDF库的方法（https://blog.csdn.net/zhebushibiaoshifu/article/details/136451536）；关于Python语言读取.nc数据的方法，大家可以参考Python读取.nc数据并提取指定时间、经纬度维度对应的变量数值（https://blog.csdn.net/zhebushibiaoshifu/article/details/136264348）。

  首先，在C++ 语言的代码读取.nc格式文件时，出现了如下图所示的情况，可以看到这些值都是负值；而实际上我这里的这个.nc格式文件肯定不应该如此。

  正常情况下，在ArcMap软件中打开上述这个.nc格式的文件，其数值正常范围的区间应该是如下图所示，肯定都是在大于0的区间内；当然，数据中确实可能会有NoData值，但尽管如此，这个.nc格式文件也不可能像上图那样，出现这么多不同的负数值。

  那么，如果出现类似上述这样的情况，大家就可以多多注意，很可能是由于存在scale和offset导致的问题了。

  首先，什么是scale和offset呢？简单来说，为了存储方便，.nc格式文件在保存数据的时候，可能会让原本的真实数据先乘以某个数，然后再加上某个数（很多.tif格式的遥感影像也是这么存储的，也就是常说的缩放系数）。例如，假设一个.nc格式文件原本的数值都是大于0、小于1的数值（例如反射率数据，都是0.X的数据），那么直接存储小数就需要占用大量的存储空间（因为需要float格式或者double格式）；而如果让这些数据都乘上1000或者10000，也就是尽可能让小数部分消除，那么就可以用int格式来存储数据，从而降低了对存储空间的占用。

  因此，如果我们待读取的.nc格式文件含有这个scale和offset，那么在使用C++ 语言中的netCDF库读取.nc格式文件时，读到的数据就是经过缩放处理后的数据；对此，我们需要手动将这个缩放后的数据，先乘上scale，再加上offset，从而得到最终的真实结果数据。这一个步骤，在Python语言的netCDF库中，应该是会自动帮我们处理（好像是这样的，因为之前用Python语言读取.nc格式文件的时候，都没有注意到过这个scale和offset）；而在C++ 语言的netCDF库中，就需要我们自行手动处理了。

  在netCDF库的官方网站中，也有关于这个scale和offset的说明——如下图所示，二者在其中分别写作scale_factor和add_offset；在官方网站中提到，只要在.nc格式文件中看到这2个参数，都需要在读取数据后，自行手动将其乘以或添加到原数据中。

  因此，在用C++ 语言netCDF库读取.nc格式的栅格文件时，如果我们是第一次读取它，那么可以通过如下的代码，获取其变量的属性。

	NcFile file(path, NcFile::read);
	NcVar var = file.getVar(type);
	map<string, NcVarAtt> attributes_map = var.getAtts();

  其中，NcFile file(path, NcFile::read);含义为创建一个NcFile对象，path是要打开的.nc格式的栅格文件的路径，NcFile::read表示以只读模式打开文件；随后，NcVar var = file.getVar(type);表示调用file对象的getVar()方法，获取了指定变量名type（也就是我们需要读取的变量）的NcVar对象；最后，map<string, NcVarAtt> attributes_map = var.getAtts();调用var对象的getAtts()方法，获取了变量的所有属性，并将它们存储在一个map<string, NcVarAtt>对象中。在这个map中，属性的名称是键，对应的NcVarAtt对象是值。

  其中，这个attributes_map如下图所示；可以看到，其中是具有scale_factor和add_offset的。

  但是，如果此时我们直接查看这个attributes_map，是看不到scale_factor和add_offset具体的值的，因为它的值还是一个NcAtt对象；如下图所示。

  我们需要通过如下的代码，首先通过.getAtt()方法获取这个属性，然后用.getValues()方法获取这个属性的具体数值。

	NcVarAtt attribute_offset = var.getAtt("add_offset");
	NcVarAtt attribute_scale = var.getAtt("scale_factor");
	double offset, scale;
	attribute_offset.getValues(&offset);
	attribute_scale.getValues(&scale);

  其中，对于上述代码，如果大家对变量值的精度有较高要求，记得要选择double类型的变量来存储scale_factor和add_offset——如果选择的是float，可能会丢失一些精度。

  运行上述代码，我们将得到如下图所示的结果。

  可以看到，scale_factor和add_offset的值都已经显示出来了。

  那么，我们就可以将这个scale_factor和add_offset，分别作用到我们读取得到的原始数据上（因为我这里.nc格式数据的数据量非常大，所以我们就只处理前100个），来看看其数值是否正确；具体代码如下。

	vector<double> var_array(time_size * latitude_size * longitude_size);
	var.getVar(var_array.data());
	for (int i = 0; i < 100; ++i) {
		var_array[i] *= scale;
		var_array[i] += offset;
	}

  可以看到，此时得到的结果，就符合实际了；如下图所示。

  此外，我们可以在ArcGIS软件中打开这个.nc格式的数据，找到其左上角的像元，获取一下这个像元的数值，如下图所示。

  可以看到，此时上图中所显示的数据，就和上上图中，我们在Visual Studio的代码中读取到的.nc文件的数据是一致的了。

  当然，这里也需要注意，有些.nc格式的数据，其变量也可能不含有scale_factor和add_offset这两个属性的，如下图所示；所以我们都可以用本文前述的代码，先获取其属性，看看到底有没有scale_factor和add_offset；如果有的话，在执行对应的数据恢复操作即可。

  至此，大功告成。

欢迎关注：疯狂学习GIS