本文介绍了如何使用netCDF4和gdal库读取和转换气象数据的NC文件为TIF格式，以及如何批量处理这些TIF文件，计算年累计蒸发量。首先，通过glob库读取NC文件，然后利用netCDF4库解析文件内容，创建NC_to_tiffs函数将数据转换为TIF。接着，对TIF文件进行裁剪操作，最后，通过读取和累加裁剪后的TIF文件数据，统计年蒸发总量。整个过程涉及到数据处理和GIS操作的关键步骤。
大部分气象数据都是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.0
tavg definition:: past 3-hour average
acc definition:: past 3-hour accumulation
inst definition:: instantaneous
title: GLDAS2.1 LIS land surface model output monthly mean
institution: NASA GSFC
source: Noah_v3.6 forced with GDAS-AGRMET-GPCPv13rA1_202101
history: created on date: 2021-04-23T11:19:32.357
references: Rodell_etal_BAMS_2004, Kumar_etal_EMS_2006, Peters-Lidard_etal_ISSE_2007
conventions: CF-1.6
comment: website: https://ldas.gsfc.nasa.gov/gldas, https://lis.gsfc.nasa.gov/
MAP_PROJECTION: EQUIDISTANT CYLINDRICAL
SOUTH_WEST_CORNER_LAT: -59.875
SOUTH_WEST_CORNER_LON: -179.875
DX: 0.25
DY: 0.25
dimensions(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 内存大小M
format='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_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
elif 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
扫一扫
215
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
