利用python如何处理nc数据详解

前言

这两天帮一个朋友处理了些 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 查看数据信息,执行:

gdalinfo name.nc

即可得到如下信息:

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 的具体信息,可以执行:

gdalinfo NETCDF:name.nc:SUBDATASET_NAME

此处的 SUBDATASET_NAME 为上面的 T2、PSFC 等等,可以得到如下信息:

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 (本地下载))。读取方式如下:

dataset = netCDF4.Dataset('name.nc') # open the dataset

这样即可读出整个 nc 中的数据信息,如果需要获取某个 SUBDATASET 只需要使用 dataset[SUBDATASET_NAME] 即可,返回的是一个三维数组,表示不同时间段(或其他区分方式下)的数据信息。

我们可以对此数组做各种操作,如求平均值、方差等等,又让我想起了大学里的那一堆枯燥但又让人很有兴趣的实验课程。当然,此处如果使用 numpy 框架进行处理,会起到事半功倍的效果,如求长时间序列下的平均值:

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 数据,如下:

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 中(直白的说就是添加空间信息)。

也很简单,如下即可:

with rio.open(newfile, 'w', **out_meta) as dst:
 dst.write_band(1, res_arr)

newfile 为存储路径,res_arr 为计算结果数组,注意尺寸不要发生变化(width*height),out_meta 为目标文件的元数据描述信息,可以直接将上面 src.meta 进行简单处理即可。

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 数据。每个目标都有多条路可以达到,重要的是找到那条自己喜欢的和适合自己的路,然而话又说回来,即使走的不是想要的那条路,不是一样可以达到目标嘛!所以关键是要找到自己的目标。

好了,以上就是这篇文章的全部内容了,希望本文的内容对大家的学习或者工作具有一定的参考学习价值,如果有疑问大家可以留言交流,谢谢大家对我们的支持。

时间: 2018-05-21

MySQL数据库设计之利用Python操作Schema方法详解

弓在箭要射出之前,低声对箭说道,"你的自由是我的".Schema如箭,弓似Python,选择Python,是Schema最大的自由.而自由应是一个能使自己变得更好的机会. Schema是什么? 不管我们做什么应用,只要和用户输入打交道,就有一个原则--永远不要相信用户的输入数据.意味着我们要对用户输入进行严格的验证,web开发时一般输入数据都以JSON形式发送到后端API,API要对输入数据做验证.一般我都是加很多判断,各种if,导致代码很丑陋,能不能有一种方式比较优雅的验证用户数据呢

利用Python破解验证码实例详解

一.前言 本实验将通过一个简单的例子来讲解破解验证码的原理,将学习和实践以下知识点: Python基本知识 PIL模块的使用 二.实例详解 安装 pillow(PIL)库: $ sudo apt-get update $ sudo apt-get install python-dev $ sudo apt-get install libtiff5-dev libjpeg8-dev zlib1g-dev \ libfreetype6-dev liblcms2-dev libwebp-dev tcl

如何利用python读取micaps文件详解

最近用编程处理文件挺多的,matlab用得比较熟,但还是想用python来写写,Fortran就不用了. 所用到的数据如下图,前面4行是说明,实际要用的数据是第5行开始. 一共是有29*53个点,每一组就有53个数据,一共是有29组. 下面就是操作了 # 导入所需的库 import numpy # 打开 micaps 文件 f1 = open('13052520.000', 'rt') f2 = open('data.txt', 'wt') # 前面4行为注释数据,没有用 for i in ra

python异步存储数据详解

在Python中,数据存储方式分为同步存储和异步存储.同步写入速度比较慢,而爬虫速度比较快,有可能导致数据保存不完整,一部分数据没有入库.而异步可以将爬虫和写入数据库操作分开执行,互不影响,所以写入速度比较快,能够保证数据的完整性. 异步存储数据库大致看分为以下步骤: 1. 在settings中配置Mysql链接需要的参数(主机地址.用户账号.密码.需要操作的表名.编码格式等) 2. 自定义Pipeline,实现from_settings函数 3. from twisted.enterprise

如何利用Python模拟GitHub登录详解

前言 最近学习了Fiddler抓包工具的简单使用,通过抓包,我们可以抓取到HTTP请求,并对其进行分析.现在我准备尝试着结合Python来模拟GitHub登录. Fiddler抓包分析 首先,我们想要模拟一个网站的登录,我们必须要简单了解其大致过程. 在这里,我通过Fiddler来抓取GitHub登录的请求,从网页上登录的URL为:https://github.com/login ,抓包结果如下: 左边的是会话列表,右边的是请求和响应的数据.一般情况下,登录都是用POST请求,因为我在左边的会话

利用Python破解斗地主残局详解

前言 相信大家都玩过斗地主,规则就不再介绍了. 直接上一张朋友圈看到的残局图: 这道题我刚看到时,曾尝试用手工来破解,每次都以为找到了农民的必胜策略时,最后都发现其实农民跑不掉.由于手工破解无法穷尽所有可能性,所以这道题究竟农民有没有妙手跑掉呢,只能通过代码来帮助我们运算了. 本文将简要讲述怎么通过代码来求解此类问题,在最后会公布残局的最后结果,并开源代码以供大家吐槽. minimax 代码的核心思想是minimax.minimax可以拆解为两部分,mini和max,分别是最小和最大的意思. 直

python实现报表自动化详解

本篇文章将介绍: xlwt 常用功能 xlrd 常用功能 xlutils 常用功能 xlwt写Excel时公式的应用 xlwt写入特定目录(路径设置) xlwt Python语言中,写入Excel文件的扩展工具.可以实现指定表单.指定单元格的写入.支持excel03版到excel2013版.使用时请确保已经安装python环境 xlrd Python语言中,读取Excel的扩展工具.可以实现指定表单.指定单元格的读取.使用时请确保已经安装python环境. NOTICE: xlwt对Excel只

利用Python代码实现数据可视化的5种方法详解

前言 数据科学家并不逊色于艺术家.他们用数据可视化的方式绘画,试图展现数据内隐藏的模式或表达对数据的见解.更有趣的是,一旦接触到任何可视化的内容.数据时,人类会有更强烈的知觉.认知和交流. 数据可视化是数据科学家工作中的重要组成部分.在项目的早期阶段,你通常会进行探索性数据分析(Exploratory Data Analysis,EDA)以获取对数据的一些理解.创建可视化方法确实有助于使事情变得更加清晰易懂,特别是对于大型.高维数据集.在项目结束时,以清晰.简洁和引人注目的方式展现最终结果是非常

python时间日期函数与利用pandas进行时间序列处理详解

python标准库包含于日期(date)和时间(time)数据的数据类型,datetime.time以及calendar模块会被经常用到. datetime以毫秒形式存储日期和时间,datetime.timedelta表示两个datetime对象之间的时间差. 下面我们先简单的了解下python日期和时间数据类型及工具 给datetime对象加上或减去一个或多个timedelta,会产生一个新的对象 from datetime import datetime from datetime impo