快捷導(dǎo)航

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

更新時(shí)間：2024年03月08日 10:09:08 作者：瘋狂學(xué)習(xí)GIS

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

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

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

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

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

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

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

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

因此，在用C++ 語(yǔ)言netCDF庫(kù)讀取.nc格式的柵格文件時(shí)，如果我們是第一次讀取它，那么可以通過如下的代碼，獲取其變量的屬性。

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

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

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

但是，如果此時(shí)我們直接查看這個(gè)attributes_map，是看不到scale_factor和add_offset具體的值的，因?yàn)樗闹颠€是一個(gè)NcAtt對(duì)象；如下圖所示。

我們需要通過如下的代碼，首先通過.getAtt()方法獲取這個(gè)屬性，然后用.getValues()方法獲取這個(gè)屬性的具體數(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);

其中，對(duì)于上述代碼，如果大家對(duì)變量值的精度有較高要求，記得要選擇double類型的變量來(lái)存儲(chǔ)scale_factor和add_offset——如果選擇的是float，可能會(huì)丟失一些精度。

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

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

那么，我們就可以將這個(gè)scale_factor和add_offset，分別作用到我們讀取得到的原始數(shù)據(jù)上（因?yàn)槲疫@里.nc格式數(shù)據(jù)的數(shù)據(jù)量非常大，所以我們就只處理前100個(gè)），來(lái)看看其數(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;
	}

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

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

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

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