气象数据NC批量转TIF及月尺度计算年尺度-同理其他指标

大部分气象数据都是NC、NC4文件格式，因此如何读取NC文件，对于气象数据而言是十分重要的。本章主要涉及NC转TIF计算，以及批量读取TIF，批量计算年累计蒸发。

涉及库包括netCDF4、gdal，安装库的流程内容在本文中不多叙述。

导入库

import glob
import netCDF4 as nc
import numpy as np
from osgeo import gdal,osr,ogr

第一步读取文件夹文件及读取nc文件

glob是一个很好用的读取文件夹中特定文件的库

Input_folder = r'.\nc'
Output_folder = r'.\tif\output'
# 读取所有nc数据
data_list = glob.glob(Input_folder + '/*.nc4')
data_list[:5]

out[2]：
['.\\nc\\GLDAS_NOAH025_M.A202101.021.nc4','.\\nc\\GLDAS_NOAH025_M.A202102.021.nc4','.\\nc\\GLDAS_NOAH025_M.A202103.021.nc4','.\\nc\\GLDAS_NOAH025_M.A202104.021.nc4','.\\nc\\GLDAS_NOAH025_M.A202105.021.nc4']

读取nc文件基本内容，关于注释的两行内容，读者可以自行查看

具体有几项需要说明一下，后续定义NC_to_tiffs函数时需要修改部分内容，包括缺失值missing_value，变量名称variables，dimensions中对于经纬度的名称定义lons，lats

nc_data_obj = nc.Dataset(data_list[0])
print(nc_data_obj, type(nc_data_obj))  # 了解NC_DS的数据类型，<class 'netCDF4._netCDF4.Dataset'>
# print(nc_data_obj.variables)  # 了解变量的基本信息
# print(nc_data_obj)

out[3]：
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):missing_value: -9999.0tavg definition:: past 3-hour averageacc definition:: past 3-hour accumulationinst definition:: instantaneoustitle: GLDAS2.1 LIS land surface model output monthly meaninstitution: NASA GSFCsource: Noah_v3.6 forced with GDAS-AGRMET-GPCPv13rA1_202101history: created on date: 2021-04-23T11:19:32.357references: Rodell_etal_BAMS_2004, Kumar_etal_EMS_2006, Peters-Lidard_etal_ISSE_2007conventions: CF-1.6comment: website: https://ldas.gsfc.nasa.gov/gldas, https://lis.gsfc.nasa.gov/MAP_PROJECTION: EQUIDISTANT CYLINDRICALSOUTH_WEST_CORNER_LAT: -59.875SOUTH_WEST_CORNER_LON: -179.875DX: 0.25DY: 0.25dimensions(sizes): lon(1440), lat(600), time(1), bnds(2)variables(dimensions): float32 lat(lat), float32 lon(lon), float64 time(time), float64 time_bnds(time, bnds), float32 Swnet_tavg(time, lat, lon), float32 Lwnet_tavg(time, lat, lon), float32 Qle_tavg(time, lat, lon), float32 Qh_tavg(time, lat, lon), float32 Qg_tavg(time, lat, lon), float32 Snowf_tavg(time, lat, lon), float32 Rainf_tavg(time, lat, lon), float32 Evap_tavg(time, lat, lon), float32 Qs_acc(time, lat, lon), float32 Qsb_acc(time, lat, lon), float32 Qsm_acc(time, lat, lon), float32 AvgSurfT_inst(time, lat, lon), float32 Albedo_inst(time, lat, lon), float32 SWE_inst(time, lat, lon), float32 SnowDepth_inst(time, lat, lon), float32 SoilMoi0_10cm_inst(time, lat, lon), float32 SoilMoi10_40cm_inst(time, lat, lon), float32 SoilMoi40_100cm_inst(time, lat, lon), float32 SoilMoi100_200cm_inst(time, lat, lon), float32 SoilTMP0_10cm_inst(time, lat, lon), float32 SoilTMP10_40cm_inst(time, lat, lon), float32 SoilTMP40_100cm_inst(time, lat, lon), float32 SoilTMP100_200cm_inst(time, lat, lon), float32 PotEvap_tavg(time, lat, lon), float32 ECanop_tavg(time, lat, lon), float32 Tveg_tavg(time, lat, lon), float32 ESoil_tavg(time, lat, lon), float32 RootMoist_inst(time, lat, lon), float32 CanopInt_inst(time, lat, lon), float32 Wind_f_inst(time, lat, lon), float32 Rainf_f_tavg(time, lat, lon), float32 Tair_f_inst(time, lat, lon), float32 Qair_f_inst(time, lat, lon), float32 Psurf_f_inst(time, lat, lon), float32 SWdown_f_tavg(time, lat, lon), float32 LWdown_f_tavg(time, lat, lon)groups:  <class 'netCDF4._netCDF4.Dataset'>

定义NC_to_tiffs函数

继续重复，对于不同nc文件需要注意修改缺失值missing_value，变量名称variables，dimensions中对于经纬度的名称定义lons，lats。

如果lats是从小到大排序，需要设置 u_arr[ i ] [::-1]

另外需要考虑数据单位问题，在本章节中的蒸发量为例，GLDAS对于 “Evap_tavg” 的单位定义为 kg m-2 s-1 ，kg m-2 相当于 mm，所以对于月尺度的GLDAS数据蒸发量而言，只需要乘以时间即可。

def NC_to_tiffs(data, Output_folder, bandname):nc_data_obj = nc.Dataset(data)
#     print(nc_data_obj, type(nc_data_obj))  # 了解NC_DS的数据类型，<class 'netCDF4._netCDF4.Dataset'>
#     print(nc_data_obj.variables)  # 了解变量的基本信息
#     print(nc_data_obj)Lon = nc_data_obj.variables['lon'][:]Lat = nc_data_obj.variables['lat'][:]u_arr = np.asarray(nc_data_obj.variables[bandname])  # 这里根据需求输入想要转换的波段名称# 影像的左上角和右下角坐标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(u_arr[:])):# 创建.tif文件driver = gdal.GetDriverByName('GTiff')out_tif_name = Output_folder + '/' + bandname +'_'  + data.split('\\')[2].split('.')[1][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())  # 给新建图层赋予投影信息# 去除异常值u_arr[u_arr[:, :] == -9999.0] = np.nan# 数据写出out_tif.GetRasterBand(1).WriteArray(u_arr[i][::-1]*60*60*24*30)# 将数据写入内存，此时没有写入硬盘 此处[::-1]用于图像的垂直镜像对称，避免图像颠倒out_tif.FlushCache()  # 将数据写入硬盘del out_tif  # 注意必须关闭tif文件

导出TIF并保存

for i in range(len(data_list)):data = data_list[i]NC_to_tiffs(data, Output_folder,"Evap_tavg")

裁剪TIF

虽然获得了全球蒸发量数据，换到GIS软件中裁剪也是可以的，但是好麻烦，尤其是遇到需要统计年单位结果时候。

读取TIF数据

Input_folder = './tif/output/'
Output_folder = './tif/clip/'
# 读取所有nc数据
data_list = glob.glob(Input_folder + '/*.tif')
data_list[:5]

定义裁剪函数

def clipTif(data):input_shape = r"shp/Uzb.shp" output_raster= Output_folder + data.split('\\')[1].split('.')[0].split("_")[-1] +".tif"# tif输入路径，打开文件input_raster = data# 矢量文件路径，打开矢量文件input_raster=gdal.Open(input_raster)ds = gdal.Warp(output_raster, # 裁剪图像保存完整路径（包括文件名）input_raster, # 待裁剪的影像# warpMemoryLimit=500 内存大小Mformat='GTiff', # 保存图像的格式cutlineDSName = input_shape, # 矢量文件的完整路径cropToCutline = True,copyMetadata=True,creationOptions=['COMPRESS=LZW',"TILED=True"],cutlineWhere="FIELD = 'whatever'",dstNodata=None)

裁剪

for i in range(len(data_list)):data = data_list[i]clipTif(data)

统计年蒸发总量

定义读取tif函数

def readTif(fileName):dataset = gdal.Open(fileName)im_width = dataset.RasterXSize #栅格矩阵的列数im_height = dataset.RasterYSize #栅格矩阵的行数
#     print(im_width, im_height)im_data = dataset.ReadAsArray(0,0,im_width,im_height)#获取数据im =  im_data[0:im_height,0:im_width]#获取蓝波段return im

读取文件夹

Input_folder = './tif/clip/'
Output_folder = './tif/sum/'
# 读取所有nc数据
data_list = glob.glob(Input_folder + '/*.tif')
data_list[:5]

设置输出的基本格式[方便批处理时候的标准参考]

dataset = gdal.Open(data_list[0])
im_width = dataset.RasterXSize #栅格矩阵的列数
im_height = dataset.RasterYSize #栅格矩阵的行数
im_bands = dataset.RasterCount #波段数
im_data = dataset.ReadAsArray(0,0,im_width,im_height)#获取数据
im_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息
im_proj = dataset.GetProjection()#获取投影信息

定义写入函数

def writeTiff(im_data,path,im_width,im_height,im_bands,im_geotrans,im_proj):if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelif len(im_data.shape) == 2:im_data = np.array([im_data])else:im_bands, (im_height, im_width) = 1,im_data.shape#创建文件driver = gdal.GetDriverByName("GTiff")path = Output_folder + "\\" + path.split('\\')[1].split('.')[0][:4] + ".tif"dataset = driver.Create(path, im_width, im_height, im_bands, datatype)if(dataset!= None):dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数dataset.SetProjection(im_proj) #写入投影for i in range(im_bands):dataset.GetRasterBand(i+1).WriteArray(im_data[i])del dataset

输出年结果

需要注意的问题是年单位时候对于降水、蒸发等的定义问题，可以对tmp进行修改

index = 1
tmp = np.zeros(readTif(data_list[0]).shape ,dtype=float)
for i in range(0, len(data_list)):if index % 12 == 0:print(index / 12)writeTiff(tmp,data_list[i -1] ,im_width ,im_height ,im_bands ,im_geotrans ,im_proj)tmp = np.zeros(readTif(data_list[i]).shape ,dtype=float)tmp += np.array(readTif(data_list[i]))index+=1