基于Python+GDAL实现nc格式转geotiff格式

1. 目的
2. 版本
3. 基础知识
- 3.1 什么是nc文件？
- 3.2 基于Python处理nc文件需要用到的库
4. 程序示例
5. 问题
- 5.1 影像分辨率的设置
参考资料

1. 目的

（1）掌握基于Python处理nc格式文件的基本方法
（2）学会将程序函数化，提高程序可读性【待补充】

2. 版本

（1）2019年6月20日; Version 1

3. 基础知识

3.1 什么是nc文件？

NetCDF(network Common Data Form)网络通用数据格式是由美国大学大气研究协会（University Corporation for Atmospheric Research，UCAR)的Unidata项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。用户可以借助多种方式方便地管理和操作NetCDF 数据集^[1]。

3.2 基于Python处理nc文件需要用到的库

基于Python处理nc数据必要的库是netCDF4，同时如果需要将nc文件转换为tiff文件，还需要osgeo库中的gdal子库和osr子库。

4. 程序示例

以Python将nc格式的标准化降水蒸散发指数（SPEI）文件转换为tif格式的文件为例，说明nc转tif的python实现过程，SPEI数据下载自西班牙国家研究委员会（CSIC）下属的SPEI下载网站^[2]。


'''1. 目的：基于Python将.nc文件转换为.tif文件2. 山东青岛  2019年6月20日
'''
# 0. 模块导入
import numpy as np
import netCDF4 as Dataset
from osgeo import gdal,osr# 1. 路径处理和变量定义
RootDir = r'F:\SPEI_NC'
SPEI_NC = RootDir + '\\spei12.nc' # 输入文件
# 1.1 输出路径-OutPath
OutPath = RootDir + '\\SPEI_TIF'
if os.path.exists(OutPath):shutil.rmtree(OutPath)os.mkdir(OutPath)
else:os.mkdir(OutPath)
OutTif = OutPath + '\\Global_SPEI_Test.tif' # 2. NetCDF文件处理
NC_DS = Dataset(SPEI_NC)
print(NC_DS,type(NC_DS)) # 了解NC_DS的数据类型，<class 'netCDF4._netCDF4.Dataset'>
print(NC_DS.variables) # 了解变量的基本信息
print(NC_DS.variables['spei']) # 了解SPEI的基本信息
Lat = NC_DS.variables['lat'][:]
Lon = NC_DS.variables['Lon'][:]
SPEI = NC_DS.variables['spei'][13,:,:] # 1901年1月的SPEI_12，注意SPEI的存储形式决定了读取方式
print(type(SPEI),SPEI.shape) # 了解SPEI的数据类型，和维数# 2.1 异常值处理
SPEI = np.asarray(SPEI) # 数据类型转换
SPEI[np.where(SPEI == 1.e+30)] = -99 # 3. 将数据写出到.tif文件中# 3.1 影像的左上角和右下角坐标
LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] # 3.2 影像的分辨率，此处float(N_Lon)-1是为了保证分辨率为0.5 degree，不知是否合理，望指正
N_Lat = len(Lat)
N_Lon = len(Lon)
Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)# 3.3 构建.tiff文件框架
spei_ds = gdal.GetDriverByName('Gtiff').Create(OutTif,N_Lon,N_Lat,1,gdal.GDT_Float32) # 3.4 设置影像的显示范围
geotransform = (LonMin,Lon_Res, 0, LatMin, 0, Lat_Res)
spei_ds.SetGeoTransform(geotransform)# 3.5 地理坐标系统信息
srs = osr.SpatialReference() #获取地理坐标系统信息，用于选取需要的地理坐标系统
print(type(srs))
print(srs)
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
spei_ds.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息# 3.6 数据写出
spei_ds.GetRasterBand(1).WriteArray(SPEI) # 将数据写入内存，此时没有写入硬盘
spei_ds.FlushCache() # 将数据写入硬盘
spei_ds = None # 关闭spei_ds指针，注意必须关闭print('Finished')

5. 问题

5.1 影像分辨率的设置

 # 3.2 影像的分辨率，此处float(N_Lon)-1是为了保证分辨率为0.5 degree，不知是否合理，望指正
N_Lat = len(Lat)
N_Lon = len(Lon)
Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)

参考资料

[1] : netCDF百度百科
[2] : MATLAB中利用ncread函数读取nc文件