Python读取.nc数据并提取指定时间、经纬度维度对应的变量数值
疯狂学习GIS 2024-07-27 08:35:09 阅读 100
本文介绍基于Python语言的<code>netCDF4库,读取.nc
格式的数据文件,并提取指定维(时间、经度与纬度)下的变量数据的方法。
我们之前介绍过.nc
格式的数据,其是NetCDF
(Network Common Data Form)文件的扩展名,是一种常用的科学数据存储格式,多用于存储科学和工程领域的大型数据集。同时,在我们之前的文章Python批量读取大量nc格式文件并导出全部时间信息(https://blog.csdn.net/zhebushibiaoshifu/article/details/135331417)中,就介绍过基于netCDF4
库,对一个文件夹下大量.nc
格式数据文件的某一维的信息加以提取的方法。而在本文中,我们则是同样基于netCDF4
库,读取.nc
文件,并提取指定维(Dimensions,也就相当于是自变量)下的变量(Variables,也就相当于是因变量)的具体数值。
首先,我们需要配置一下netCDF4
库,具体配置方法大家可以参考文章Anaconda下Python中h5py与netCDF4模块下载与安装方法(https://blog.csdn.net/zhebushibiaoshifu/article/details/120553597)。
随后,本文所需代码如下。
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 22 21:41:52 2024
@author: fkxxgis
"""
import netCDF4 as nc
nc_path = r"F:\Data_Reflectance_Rec\soil_1\2020_01.nc"
nc_data = nc.Dataset(nc_path)
print(nc_data)
time_value = nc_data.variables["time"][:]
longitude_value = nc_data.variables["lon"][:]
latitude_value = nc_data.variables["lat"][:]
# 第一种需求
time_need = 0
nc_value_1 = nc_data.variables["swvl1"][time_need, : , : ]
# 第二种需求
longitude_need = 106.467
latitude_need = 36.817
longitude_nc = (abs(longitude_value - longitude_need)).argmin()
latitude_nc = (abs(latitude_value - latitude_need)).argmin()
nc_value_2 = nc_data.variables["swvl1"][time_need, latitude_nc, longitude_nc]
其中,我们首先导入netCDF4
库,并指定要读取的.nc
格式数据文件的路径nc_path
;随后,使用nc.Dataset()
打开这一文件,并将返回的Dataset
对象存储在nc_data
变量中;紧接着,通过print()
打印nc_data
,这将显示要读取的.nc
格式数据文件的基本信息,如变量、维、属性等——这里具体打印出来的情况如下图所示。
其中,在上图我们需要重点关注紫色框内的内容。首先,在<code>dimensions中,我们可以看到所有的维;我这里的.nc
格式数据是一个表示气象的数据,所以文件中的维依次就是时间、纬度与经度;随后,在variables
中,我们可以看到所有的数据变量(这里的数据变量是包含了维、变量与其他参数)——其中我们重点观察数据中的因变量(也就是上图中的swvl1
),需要留意一下其后不同维的排序顺序,在后面我们按照维提取变量数据的时候会用到。
回到前述代码的介绍中。通过前面print()
打印出来的nc_data
信息,我们知道了这个.nc
数据的维,此时我们可以将这几个维也打印出来看看。例如,time_value = nc_data.variables["time"][:]
就表示将时间这个维打印出来,相当于获取了全部的时间节点。
再接下来,我们即可开始按照维来提取变量。为了方便,我们就以这个.nc
文件的时间维中的第一个节点对应的数据(也就是第一景数据)为例来介绍;因此,我们先将time_need
设置为0
,表示读取第一个时间节点的数据。在这里,我们给出了2
种按照维来提取变量的需求。
首先,是第一种需求,也就是time_need = 0
这一行代码的下一行。nc_data.variables["swvl1"]
表示这个.nc
文件中读取名为swvl1
的变量的值,而后面的[time_need, :, :]
表示选择指定时间维下的所有经度和纬度位置的值。这些值将被存储在nc_value_1
变量中,也就是说这个nc_value_1
变量相当于就是当前这个.nc
文件的第一景数据(时间节点排在第一位的数据)。
其次,是第二种需求。前面我们提取了指定时间维下的所有经度和纬度位置的值,那么现在就更进一步,提取指定时间维度、经度维度以及纬度维度的数据(相当于就是从前面的一景数据变成了一个像元的数据)。首先,我们指定一个处于.nc
文件成像范围内的目标经度longitude_need
和目标纬度latitude_need
,并使用argmin()
函数找到目标经、纬度值与文件中经度、纬度的维数据值最接近的索引值——即longitude_nc
和latitude_nc
。最后,即可使用nc_data.variables["swvl1"][time_need, latitude_nc, longitude_nc]
来获取特定时间、经度和纬度位置的值,并将结果存储在nc_value_2
变量中。
这里提一句——为什么需要用argmin()
函数呢?这个是因为,我们在实际情况中,需要提取指定空间位置的像素时,这个位置的经、纬度数据肯定是随机的;而通过argmin()
函数,就可以找到.nc
文件里面经度、纬度所对应的维中,与我们实际需要的经、纬度最接近的那个数值所对应的维的下标。例如,上述代码中,我们希望提取实际经度为106.467
位置处的数据;而我这里这个.nc
文件,其维中的经度的分辨率是0.1
,那它自然没有办法非常精确地确定106.467
的位置;所以需要通过argmin()
函数,找到与106.467
最接近的数据106.5
,并进一步确定出这个106.5
所在的经度维的下标,那么就可以提取出指定的变量了。
如下图所示,我们通过上述第二种需求,提取出来了目标时间、经度与纬度维下的一个像素。
那么这个像素值对不对呢?我们可以在ArcGIS中打开这个<code>.nc文件的第一景数据,找到代码中目标经、纬度(也就是longitude_need
和latitude_need
所指向的数据)所对应的像元,并查看其像素值;如下图所示。
可以看到,上图在ArcGIS提取出来的像素值,与上上图在Python中提取出来的像素值一致,说明我们的代码无误。
至此,大功告成。
欢迎关注:疯狂学习GIS
声明
本文内容仅代表作者观点,或转载于其他网站,本站不以此文作为商业用途
如有涉及侵权,请联系本站进行删除
转载本站原创文章,请注明来源及作者。