C++讀取NC數(shù)據(jù)的結(jié)果與真實數(shù)值不一致的解決方法

更新時間：2024年03月08日 10:09:08 作者：瘋狂學習GIS

本文介紹基于C++ 語言的netCDF庫讀取.nc格式的柵格文件時,代碼讀取到的數(shù)據(jù)與柵格文件的實際數(shù)據(jù)不一致的解決方法,文中通過代碼示例和圖文講解的非常詳細,需要的朋友可以參考下

最近，由于需要讀取ERA5氣象數(shù)據(jù)，因此使用C++ 語言中的netCDF庫讀取.nc格式文件。其中，偶然發(fā)現(xiàn)在Visual Studio的代碼中讀取到的.nc文件的數(shù)據(jù)，和其實際的數(shù)據(jù)（ArcMap等軟件打開所顯示的數(shù)據(jù)）不一致；這里就介紹一種可能導致上述情況的原因，以及對應的解決方法。

首先，在C++ 語言的代碼讀取.nc格式文件時，出現(xiàn)了如下圖所示的情況，可以看到這些值都是負值；而實際上我這里的這個.nc格式文件肯定不應該如此。

正常情況下，在ArcMap軟件中打開上述這個.nc格式的文件，其數(shù)值正常范圍的區(qū)間應該是如下圖所示，肯定都是在大于0的區(qū)間內(nèi)；當然，數(shù)據(jù)中確實可能會有NoData值，但盡管如此，這個.nc格式文件也不可能像上圖那樣，出現(xiàn)這么多不同的負數(shù)值。

那么，如果出現(xiàn)類似上述這樣的情況，大家就可以多多注意，很可能是由于存在scale和offset導致的問題了。

首先，什么是scale和offset呢？簡單來說，為了存儲方便，.nc格式文件在保存數(shù)據(jù)的時候，可能會讓原本的真實數(shù)據(jù)先乘以某個數(shù)，然后再加上某個數(shù)（很多.tif格式的遙感影像也是這么存儲的，也就是常說的縮放系數(shù)）。例如，假設(shè)一個.nc格式文件原本的數(shù)值都是大于0、小于1的數(shù)值（例如反射率數(shù)據(jù)，都是0.X的數(shù)據(jù)），那么直接存儲小數(shù)就需要占用大量的存儲空間（因為需要float格式或者double格式）；而如果讓這些數(shù)據(jù)都乘上1000或者10000，也就是盡可能讓小數(shù)部分消除，那么就可以用int格式來存儲數(shù)據(jù)，從而降低了對存儲空間的占用。

因此，如果我們待讀取的.nc格式文件含有這個scale和offset，那么在使用C++ 語言中的netCDF庫讀取.nc格式文件時，讀到的數(shù)據(jù)就是經(jīng)過縮放處理后的數(shù)據(jù)；對此，我們需要手動將這個縮放后的數(shù)據(jù)，先乘上scale，再加上offset，從而得到最終的真實結(jié)果數(shù)據(jù)。這一個步驟，在Python語言的netCDF庫中，應該是會自動幫我們處理（好像是這樣的，因為之前用Python語言讀取.nc格式文件的時候，都沒有注意到過這個scale和offset）；而在C++ 語言的netCDF庫中，就需要我們自行手動處理了。

在netCDF庫的官方網(wǎng)站中，也有關(guān)于這個scale和offset的說明——如下圖所示，二者在其中分別寫作scale_factor和add_offset；在官方網(wǎng)站中提到，只要在.nc格式文件中看到這2個參數(shù)，都需要在讀取數(shù)據(jù)后，自行手動將其乘以或添加到原數(shù)據(jù)中。

因此，在用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);含義為創(chuàng)建一個NcFile對象，path是要打開的.nc格式的柵格文件的路徑，NcFile::read表示以只讀模式打開文件；隨后，NcVar var = file.getVar(type);表示調(diào)用file對象的getVar()方法，獲取了指定變量名type（也就是我們需要讀取的變量）的NcVar對象；最后，map<string, NcVarAtt> attributes_map = var.getAtts();調(diào)用var對象的getAtts()方法，獲取了變量的所有屬性，并將它們存儲在一個map<string, NcVarAtt>對象中。在這個map中，屬性的名稱是鍵，對應的NcVarAtt對象是值。

其中，這個attributes_map如下圖所示；可以看到，其中是具有scale_factor和add_offset的。

但是，如果此時我們直接查看這個attributes_map，是看不到scale_factor和add_offset具體的值的，因為它的值還是一個NcAtt對象；如下圖所示。

我們需要通過如下的代碼，首先通過.getAtt()方法獲取這個屬性，然后用.getValues()方法獲取這個屬性的具體數(shù)值。

	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，可能會丟失一些精度。

運行上述代碼，我們將得到如下圖所示的結(jié)果。

可以看到，scale_factor和add_offset的值都已經(jīng)顯示出來了。

那么，我們就可以將這個scale_factor和add_offset，分別作用到我們讀取得到的原始數(shù)據(jù)上（因為我這里.nc格式數(shù)據(jù)的數(shù)據(jù)量非常大，所以我們就只處理前100個），來看看其數(shù)值是否正確；具體代碼如下。

	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;
	}

可以看到，此時得到的結(jié)果，就符合實際了；如下圖所示。

此外，我們可以在ArcGIS軟件中打開這個.nc格式的數(shù)據(jù)，找到其左上角的像元，獲取一下這個像元的數(shù)值，如下圖所示。

可以看到，此時上圖中所顯示的數(shù)據(jù)，就和上上圖中，我們在Visual Studio的代碼中讀取到的.nc文件的數(shù)據(jù)是一致的了。

當然，這里也需要注意，有些.nc格式的數(shù)據(jù)，其變量也可能不含有scale_factor和add_offset這兩個屬性的，如下圖所示；所以我們都可以用本文前述的代碼，先獲取其屬性，看看到底有沒有scale_factor和add_offset；如果有的話，在執(zhí)行對應的數(shù)據(jù)恢復操作即可。