前言
这两天帮一个朋友处理了些 nc 数据,本以为很简单的事情,没想到里面涉及到了很多的细节和坑,无论是“知难行易”还是“知易行难”都不能充分的说明问题,还是“知行合一”来的更靠谱些,既要知道理论又要知道如何实现,于是经过不太充分的研究后总结成此文,以记录如何使用 python 处理 nc 数据。
一、nc 数据介绍
nc 全称 netCDF(The Network Common Data Form),可以用来存储一系列的数组,就是这么简单(参考https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_introduction.html)。
既然 nc 可以用来一系列的数组,所以经常被用来存储科学观测数据,最好还是长时间序列的。
试想一下一个科学家每隔一分钟采集一次实验数据并存储了下来,如果不用这种格式存储,时间长了可能就需要创建一系列的 csv 或者 txt 等,而采用 nc 一个文件就可以搞定,是不是很方便。
更方便的是如果这个科学实验与气象、水文、温度等地理信息稍微沾点边的,完全也可以用 nc 进行存储, GeoTiff 顶多能多存几个波段(此处波段可以认为是气象、水文等不同信号),而 nc 可以存储不同波段的长时间观测结果,是不是非常方便。
可以使用 gdal 查看数据信息,执行:
1
|
gdalinfo name.nc
|
即可得到如下信息:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
|
Driver: netCDF / Network Common Data Format
Files: test.nc
Size is 512 , 512
Coordinate System is `'
Subdatasets:
SUBDATASET_1_NAME = NETCDF: "test.nc" :T2
SUBDATASET_1_DESC = [ 696x130x120 ] T2 ( 32 - bit floating - point)
SUBDATASET_2_NAME = NETCDF: "test.nc" :PSFC
SUBDATASET_2_DESC = [ 696x130x120 ] PSFC ( 32 - bit floating - point)
SUBDATASET_3_NAME = NETCDF: "test.nc" :Q2
SUBDATASET_3_DESC = [ 696x130x120 ] Q2 ( 32 - bit floating - point)
SUBDATASET_4_NAME = NETCDF: "test.nc" :U10
SUBDATASET_4_DESC = [ 696x130x120 ] U10 ( 32 - bit floating - point)
SUBDATASET_5_NAME = NETCDF: "test.nc" :V10
SUBDATASET_5_DESC = [ 696x130x120 ] V10 ( 32 - bit floating - point)
SUBDATASET_6_NAME = NETCDF: "test.nc" :RAINC
SUBDATASET_6_DESC = [ 696x130x120 ] RAINC ( 32 - bit floating - point)
SUBDATASET_7_NAME = NETCDF: "test.nc" :SWDOWN
SUBDATASET_7_DESC = [ 696x130x120 ] SWDOWN ( 32 - bit floating - point)
SUBDATASET_8_NAME = NETCDF: "test.nc" :GLW
SUBDATASET_8_DESC = [ 696x130x120 ] GLW ( 32 - bit floating - point)
SUBDATASET_9_NAME = NETCDF: "test.nc" :LAT
SUBDATASET_9_DESC = [ 130x120 ] LAT ( 32 - bit floating - point)
SUBDATASET_10_NAME = NETCDF: "test.nc" : LONG
SUBDATASET_10_DESC = [ 130x120 ] LONG ( 32 - bit floating - point)
Corner Coordinates:
Upper Left ( 0.0 , 0.0 )
Lower Left ( 0.0 , 512.0 )
Upper Right ( 512.0 , 0.0 )
Lower Right ( 512.0 , 512.0 )
Center ( 256.0 , 256.0 )
|
每一个 SUBDATASET 表示记录的是一种格式的数据(气象、水文等等),如果要想查看此 SUBDATASET 的具体信息,可以执行:
1
|
gdalinfo NETCDF:name.nc:SUBDATASET_NAME
|
此处的 SUBDATASET_NAME 为上面的 T2、PSFC 等等,可以得到如下信息:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
|
Driver: netCDF / Network Common Data Format
Files: test.nc
Size is 120 , 130
Coordinate System is `'
Metadata:
LAT #description=LATITUDE, SOUTH IS NEGATIVE
LAT #FieldType=104
LAT #MemoryOrder=XY
LAT #stagger=
LAT #units=degree_north
Corner Coordinates:
Upper Left ( 0.0 , 0.0 )
Lower Left ( 0.0 , 130.0 )
Upper Right ( 120.0 , 0.0 )
Lower Right ( 120.0 , 130.0 )
Center ( 60.0 , 65.0 )
Band 1 Block = 120x1 Type = Float32, ColorInterp = Undefined
NoData Value = 9.96920996838686905e + 36
Unit Type : degree_north
Metadata:
description = LATITUDE, SOUTH IS NEGATIVE
FieldType = 104
MemoryOrder = XY
NETCDF_VARNAME = LAT
stagger =
units = degree_north
|
此处只有一个 Band ,每一个 Band 记录了一个时间点(或者其他区分形式)的一条记录,这个记录是一个数组。
所以看到这里,各位应该已经明白了,可以直接使用 GDAL 处理 nc 数据,比如直接使用 gdalwarp 将某个 SUBDATASET 转成 GeoTiff 等等,此处暂且不表,各位只需要查阅一下 gdalwarp 手册即可知道如何处理。
明白了以上信息基本也就清楚了如何处理此数据。
二、数据处理
python 是运用非常广泛,自然其下各种类库非常丰富,专业一点的说法就叫生态丰富。
2.1 netCDF4
此框架可以直接将 nc 读取成数组(详细信息参考https://github.com/Unidata/netcdf4-python )。读取方式如下:
1
|
dataset = netCDF4.Dataset( 'name.nc' ) # open the dataset
|
这样即可读出整个 nc 中的数据信息,如果需要获取某个 SUBDATASET 只需要使用 dataset[SUBDATASET_NAME] 即可,返回的是一个三维数组,表示不同时间段(或其他区分方式下)的数据信息。
我们可以对此数组做各种操作,如求平均值、方差等等,又让我想起了大学里的那一堆枯燥但又让人很有兴趣的实验课程。当然,此处如果使用 numpy 框架进行处理,会起到事半功倍的效果,如求长时间序列下的平均值:
1
2
|
np_arr = np.asarray(dataset[SUBDATASET_NAME])
average_arr = np.average(np_arr, axis = 0 )
|
到这里跟地信有关的同志都会看出一个问题,此框架只能对数据进行处理,而不能进行与位置有关的操作,这就导致数据无法变成直白的地图可视化效果。其实任何数据都是相通的,我们可以采用此种方式处理完后转为 GeoTiff 等,当然我们也可以直接采用 GeoTiff 的处理流程来进行处理。
2.2 rasterio
rasterio 是 Mapbox 开源的空间数据处理框架,功能非常强大,此处不细说,只表如何处理我们的 nc 数据。
当然第一种方式就是使用 netCDF4 处理完之后,使用此框架写入 GeoTiff,但是这样不太优雅,而且使用了两个框架,明显过于麻烦,我们直接使用此框架从读数据开始处理。
此处读的时候就有技巧了,要像采用 gdalinfo 读取 SUBDATASET 一样来直接读取此 SUBDATASET 数据,如下:
1
2
3
4
|
with rio. open ( 'NETCDF:name.nc:SUBDATASET_NAME' ) as src:
print (src.meta)
dim = int (src.meta[ 'count' ])
src.read( range ( 1 , dim + 1 ))
|
即给 open 函数传入 NETCDF:name.nc:SUBDATASET_NAME,采用 src.read(range(1, dim + 1)) 可以直接读出此范围内所有 Band (时间点)的信息,范围可以自己设定,注意从 0 开始,当然也可以仅读取某个 Band 的信息。
src.meta 记录了此 SUBDATASET 的元数据信息,与 gdalinfo 看到的基本相同。
这样我们就可以继续将此数据使用 numpy 等框架进行处理,处理完之后更重要的是要写入 GeoTiff 中(直白的说就是添加空间信息)。
也很简单,如下即可:
1
2
|
with rio. open (newfile, 'w' , * * out_meta) as dst:
dst.write_band( 1 , res_arr)
|
newfile 为存储路径,res_arr 为计算结果数组,注意尺寸不要发生变化(width*height),out_meta 为目标文件的元数据描述信息,可以直接将上面 src.meta 进行简单处理即可。
1
2
3
4
5
6
7
|
out_meta =
meta.update({ "driver" : "GTiff" ,
"dtype" : "float32" ,
'count' : 1 ,
'crs' : 'Proj4: +proj=longlat +datum=WGS84 +no_defs' ,
'transform' : rasterio.transform.from_bounds(west, south, east, north, width, height)
})
|
crs 表示目标数据空间投影信息,transform 表示目标文件 空间范围信息,可以通过经纬度信息和图像尺寸等计算得到。
dst.write_band 将数据写入对应波段,当然此处也可以写入多个波段,根据计算结果而定,同样从 1 开始。
三、总结
本文简单介绍了 nc 数据的特点及如何使用 python 处理 nc 数据。每个目标都有多条路可以达到,重要的是找到那条自己喜欢的和适合自己的路,然而话又说回来,即使走的不是想要的那条路,不是一样可以达到目标嘛!所以关键是要找到自己的目标。
好了,以上就是这篇文章的全部内容了,希望本文的内容对大家的学习或者工作具有一定的参考学习价值,如果有疑问大家可以留言交流,谢谢大家对服务器之家的支持。
原文链接:https://www.cnblogs.com/shoufengwei/p/9068379.html