先看函数
1.提取NC数据和地理信息
def nc2tiff(nc_file, output_dir):
data_name = os.path.basename(nc_file) #文件名
dataset = nc.Dataset(nc_file)
#############获取经纬度和时间变量#############
Lat = dataset.variables["latitude"][:] #latitude
Lon = dataset.variables["longitude"][:] #longitude
time = dataset.variables["time"][:] #era ei temp
##############设置地理参考##############
LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]##空间范围
Lon_Res = (LonMax - LonMin) / (float(cols) - 1) #lon分辨率计算
Lat_Res = (LatMax - LatMin) / (float(rows) - 1) #lat分辨率计算
geo = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res) #geo = (-9036842.762, 25067.525, 0, 9036842.763000002, 0, -25067.525) #globsnow # -Lat_Res一定要是-的
src_srs = osr.SpatialReference() #构造projection
src_srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84",AUTHORITY["EPSG","4326"]
src_srs_wkt = src_srs.ExportToWkt() # 给新建图层赋予投影信息
################提取变量####################
dayarray = np.array(dataset.variables['SWE_tavg'][0, :, :])
out_file = os.path.join(output_dir, "gldas" + nc_file[38:46] + ".tif")
print(out_file)
write_tiff(out_file, geo, src_srs_wkt, rows, cols, dayarray)
注:这里针对的是每日NC数据,如果一个NC数据包含多维时间数据,需要根据TIME循环
2.TIFF生成
def write_tiff(tiff_file, geo, proj, rows, cols, data_array):
driver = gdal.GetDriverByName("Gtiff")
outdataset = driver.Create(tiff_file, cols, rows, 1, gdal.GDT_Float32)
outdataset.SetGeoTransform(geo)
outdataset.SetProjection(proj)
band = outdataset.GetRasterBand(1)
band.WriteArray(data_array)
3.输入路径,运行
if __name__ == '__main__':
# # 结果存储路径
tiff_dir = r"D:\40\ei\tiff"
# nc数据路径
nc_file = r"D:\40\ei\nc2\.*tif"
ncFileList = getFileName(nc_file,2016,2019) ## 筛选特定年份的数据
# nc数据转换为tif数据
for ncFile in ncFileList:
print(ncFile)
nc2tiff(ncFile, tiff_dir)