python批量读取nc气象数据并转为tif

python批量nc数据转tif

各类地理数据中，NC格式是很常见的，然而这种格式ArcGIS是无法打开的。一旦下载的话nc也是多时序的（多年、多月等等），让我们看看如何批量操作吧。

直接上代码：

import numpy as np
import netCDF4 as nc
from osgeo import gdal,osr,ogr
import os
import glob
def NC_to_tiffs(data,Output_folder):
    nc_data_obj = nc.Dataset(data)
    Lon = nc_data_obj.variables['lon'][:]
    Lat = nc_data_obj.variables['lat'][:]
    PM = nc_data_obj.variables['di_01']
    tp_arr = np.asarray(nc_data_obj.variables['di_01'])#将tp数据读取为数组
    tp_arr_float = tp_arr.astype(float)
    LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()]
    N_Lat = len(Lat)
    N_Lon = len(Lon)
    Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)
    for i in range(len(tp_arr[:])):#创建.tif文件
        driver = gdal.GetDriverByName('GTiff')
        out_tif_name = Output_folder + '\\'+ data.split('\\')[-1].split('.')[0] + '_' + str(i+1) + '.tif'
        out_tif = driver.Create(out_tif_name,N_Lon,N_Lat,1,gdal.GDT_Float32)# 设置影像的显示范围#-Lat_Res一定要是-的
        geotransform = (LonMin,Lon_Res, 0, LatMax, 0, -Lat_Res)
        out_tif.SetGeoTransform(geotransform)#获取地理坐标系统信息，用于选取需要的地理坐标系统
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)# 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
        out_tif.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息#数据写出
        out_tif.GetRasterBand(1).WriteArray(tp_arr_float[i]) # 将数据写入内存，此时没有写入硬盘
        out_tif.FlushCache() # 将数据写入硬盘
        out_tif = None # 注意必须关闭tif文件

def main():
    Input_folder = 'D:/OneDrive/data/GPCC_AI/inputNC'
    Output_folder = 'D:/OneDrive/data/GPCC_AI/output'# 读取所有nc数据
    data_list = glob.glob(Input_folder + '\\*.nc')
    for i in range(len(data_list)):
        data = data_list[i]
        NC_to_tiffs(data,Output_folder)
        print(str(i) + ': ' + data)
main()
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
32
33
34
35
36
37
38
39

运行之后会打印转换成功的数据和序号：

0: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201301.nc
1: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201302.nc
2: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201303.nc
3: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201304.nc
4: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201305.nc
5: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201306.nc
6: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201307.nc
7: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201308.nc
8: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201309.nc
9: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201310.nc
10: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201311.nc
11: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201312.nc
12: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201401.nc
13: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201402.nc
14: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201403.nc
15: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201404.nc
16: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201405.nc
17: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201406.nc
18: D:/OneDrive/data/GPCC_AI/inputNC\GPCC_DI_201407.nc
...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

打开文件夹，TIF也转换好了，此时发现我们的图像是上下左右颠倒的？

由于不同NC维度不同，数据格式、类型、存储方式都可能不同，因此一定要真正理解代码，才能对其进行修改，使之适合自己的数据。

在本例中，数据上下左右发生了颠倒，我们只需要用数组反转函数即可：

numpy.flipud() # 上下翻转
numpy.flip() # 左右翻转
1
2

那么把代码加到哪里呢？让我们来看看函数是如何工作的：

主函数的前几行是比较关键的，中括号的内部是变量名：

Lon = nc_data_obj.variables['lon'][:]
Lat = nc_data_obj.variables['lat'][:]
PM = nc_data_obj.variables['di_01']
1
2
3

那么如何查看变量名呢？这就需要数据预读取，代码如下：

data = 'D:/OneDrive/data/GPCC_AI/GPCC_DI_201301.nc'
nc_data_obj = nc.Dataset(data)
nc_data_obj.variables
1
2
3

一般nc数据中lon和lat是必有的，只不过名字也许不同，可能有的是“longtitude”和“latitude”，我们要的数据是“di_01”。

PM = nc_data_obj.variables['di_01']
tp_arr = np.asarray(nc_data_obj.variables['di_01'])#将tp数据读取为数组
1
2

所以修改了上述两行

该函数主要是将nc读成多维数组，再循环写成tif，因此我们在数组这一步进行翻转：

for i in range(len(tp_arr[:])):#创建.tif文件
    tp_arr_float[i] = np.flipud(tp_arr_float[i])  
    driver = gdal.GetDriverByName('GTiff')
1
2
3

从预览图来看，结果正确！

接下来是比较重要的投影信息：

EPSG:4326是WGS84的地理坐标

srs.ImportFromEPSG(4326)
1

最后是输入路径和输出路径，输出路径是一空文件夹，输入路径记录了全部的nc文件

Input_folder = 'D:/OneDrive/data/GPCC_AI/inputNC'
Output_folder = 'D:/OneDrive/data/GPCC_AI/output'# 读取所有nc数据
1
2

径记录了全部的nc文件

Input_folder = 'D:/OneDrive/data/GPCC_AI/inputNC'
Output_folder = 'D:/OneDrive/data/GPCC_AI/output'# 读取所有nc数据
1
2

这样就能实现批量操作了，读成数组之后，可以进行你要做的任何操作！