一个月前写了一篇“使用Arcgis处理CMIP6(netCDF)数据的流程”后，一直收到朋友们的私信，也几乎天天都有交流。说来惭愧，后来我放弃了使用Arcgis进行处理这一条路，而是直接采用Python提取栅格点的时间序列数据，进行处理。

目前忙着降尺度的工作，可能会减少在知乎上的交流，因此在这里我会把目前掌握的部分数据提取方法分享出来，希望能帮到大家。

不用写代码的方法：

这个方法是我后面才发现的，优点就四个字：简单粗暴！,缺点是输出格式死板，需要就格式问题调整很久。

首先，你需要在电脑上装有java的环境(别皱眉，只要有环境就行，不需要coding);

然后，点击下面链接，下载panoply(这个软件也有制图的功能，并且成图操作非常简单友好)；

https://www.giss.nasa.gov/tools/panoply/​www.giss.nasa.gov

下载好之后打开这个软件，选择nc文件，右边可以看到数据描述性信息，选中要输出的变量，如PR，右键选择export csv，然后等着就行。

不知道为什么一直上传图片失败，图先欠着，下同...

输出之后的csv文件普遍非常大，我示例中最大的一个达到了6g以上。它的数据输出后是这样的，每个时间点的全球数据作为一个矩阵，每一列是同一经度下，每一行是同一纬度下，以分辨率360*720为例，这个矩阵的大小就是（360，720），它是月尺度的，从1850-2014一共165年，1980个月，所以一共有1980个这样的矩阵，综上，PR的数据形态就是（1980，360，720）。

但是这样的数据形式我们不好应用，下面介绍使用python的方法提取。

默认大家都拥有配置好的python3环境，我们需要的库有：netCDF4，numpy，pandas。

其中netCDF4用来查看数据，了解数据信息，制定提取策略。numpy和pandas用于提取数据，并按我们需要的格式输出。

主要是查看数据的shape：

for 

然后你会看到我在前面提过的(t,x,y)这个形式，可以理解为一个三维矩阵（t,(x,y)）,(x,y)代表数据的分辨率，如(360，720),就意味着一共有360*760个栅格点，对应到纬度和经度后，我们就可以知道这个数据是0.5度*0.5度的栅格，即50km分辨率。

那么现在我们要做的就是提取我们想要的数据了。很多朋友都在问我怎么按区域提取数据，我没有做这个工作，我是按照手里有的实测站点位置来检索栅格点并下载数据的，按区域下载的思路我放在最后说。

pr

首先定义一个函数用来找到我们所需经纬度点最近的栅格点：

import 

现在写一个for循环，提取出某个点的时间序列数据到一个list中:

for 

同理，通过循环嵌套就可以提取出所有的点的时间序列数据。

这时，有的朋友可能会发现输出的结果里没有时间这一列。这是因为我发现nc文件中的时间格式如果直接导出的话看起来类似于时间戳，但是转化后会发现对不上，所以我采取了自己生成日期序列的笨办法(月尺度的还好，日尺度的很麻烦，因为大多数模型都没有闰年这个概念，即没有2月29日，生成日尺度日期序列时要避开这一坑)，把list转为数据框后，加上日期序列这一index，再输出为csv即可。不过，如果做多模型，并且后续处理分析也在python中完成的话，我觉得完全没必要将中间数据输出。

如果要做月变化的对比，可以采取切片的方式，将list按12步长切片，每一片都是一年的12条数据。

def 

数据提取暂且说这么多，开头提到的按区域提取数据的思路，我的想法是：通过目标区域的外多边形坐标经纬度点集，然后做一个遍历判断，把在多边形内的点筛选出来，得到符合要求的点的索引号，然后把这些点的数据输出。但是因为这段时间太忙，还没有动手尝试过，大家有思路的话，也欢迎在评论区交流。