龙空技术网

C++打开NC栅格文件后发现数据异常的解决思路

疯狂学习GIS 15

前言:

当前大家对“文件路径异常”大概比较关切,大家都想要学习一些“文件路径异常”的相关内容。那么小编同时在网上网罗了一些关于“文件路径异常””的相关内容,希望大家能喜欢,兄弟们快快来学习一下吧!

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

  最近,由于需要读取ERA5气象数据,因此使用C++ 语言中的netCDF库读取.nc格式文件。其中,偶然发现在Visual Studio的代码中读取到的.nc文件的数据,和其实际的数据(ArcMap等软件打开所显示的数据)不一致;这里就介绍一种可能导致上述情况的原因,以及对应的解决方法。

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

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

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

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

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

  在netCDF库的官方网站中,也有关于这个scaleoffset的说明——如下图所示,二者在其中分别写作scale_factoradd_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_factoradd_offset的。

  但是,如果此时我们直接查看这个attributes_map,是看不到scale_factoradd_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_factoradd_offset——如果选择的是float,可能会丢失一些精度。

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

  可以看到,scale_factoradd_offset的值都已经显示出来了。

  那么,我们就可以将这个scale_factoradd_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_factoradd_offset这两个属性的,如下图所示;所以我们都可以用本文前述的代码,先获取其属性,看看到底有没有scale_factoradd_offset;如果有的话,在执行对应的数据恢复操作即可。

  至此,大功告成。

欢迎关注:疯狂学习GIS

标签: #文件路径异常