用python将nc转tif,可以运行但是没有结果,想问问是哪里出了问题
import os
import netCDF4 as nc
import numpy as np
from osgeo import gdal,osr,ogr
import glob
def nc2tif(data,Output_folder):
tmp_data = nc.Dataset(data) #利用.Dataset()方法读取nc数据
# print(tmp_data)
# <class 'netCDF4._netCDF4.Dataset'>
# root group (NETCDF3_CLASSIC data model, file format NETCDF3):
# dimensions(sizes): lon(7849), lat(5146), time(12)
# variables(dimensions): float64 lon(lon), float64 lat(lat), float64 time(time), int16 tmp(time, lat, lon)
# print(tmp_data.variables) #lon, lat, time, tmp
# {'lon': <class 'netCDF4._netCDF4.Variable'>
# float64 lon(lon)
# long_name: longitude
# unit: degree
# unlimited dimensions:
# current shape = (7849,)
# filling on, default _FillValue of 9.969209968386869e+36 used, 'lat': <class 'netCDF4._netCDF4.Variable'>
# float64 lat(lat)
# long_name: latitude
# unit: degree
# unlimited dimensions:
# current shape = (5146,)
# filling on, default _FillValue of 9.969209968386869e+36 used, 'time': <class 'netCDF4._netCDF4.Variable'>
# float64 time(time)
# long_name: time
# unit: data.01-data.12
# calendar: gregorian
# unlimited dimensions:
# current shape = (12,)
# filling on, default _FillValue of 9.969209968386869e+36 used, 'tmp': <class 'netCDF4._netCDF4.Variable'>
# int16 tmp(time, lat, lon)
# long_name: 0.1 monthly mean temperature
# unit: degree centigrade
# missing_value: -32768.0
# unlimited dimensions:
# current shape = (12, 5146, 7849)
# filling on, default _FillValue of -32767 used}
Lat_data = tmp_data.variables['lat'][:]
Lon_data = tmp_data.variables['lon'][:]
# print(Lat_data)
# [58.63129069 58.62295736 58.61462403 ... 15.77295736 15.76462403
# 15.75629069]
# print(Lon_data)
# [ 71.29005534 71.29838867 71.30672201 ... 136.67338867 136.68172201
# 136.69005534]
tmp_arr = np.asarray(tmp_data.variables['tmp'])
#影像的左上角&右下角坐标
Lonmin, Latmax, Lonmax, Latmin = [Lon_data.min(), Lat_data.max(), Lon_data.max(), Lat_data.min()]
# Lonmin, Latmax, Lonmax, Latmin
# (71.29005533854166, 58.63129069182766, 136.6900553385789, 15.756290691830095)
#分辨率计算
Num_lat = len(Lat_data) #5146
Num_lon = len(Lon_data) #7849
Lat_res = (Latmax - Latmin) / (float(Num_lat) - 1)
Lon_res = (Lonmax - Lonmin) / (float(Num_lon) - 1)
#print(Num_lat, Num_lon)
#print(Lat_res, Lon_res)
# 5146 7849
# 0.00833333333333286 0.008333333333338078
for i in range(len(tmp_arr[:])):
#i=0,1,2,3,4,5,6,7,8,9,...
#创建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, Num_lon, Num_lat, 1, gdal.GDT_Int16)
#设置影像的显示范围
#Lat_re前需要添加负号
geotransform = (Lonmin, Lon_res, 0.0, Latmax, 0.0, -Lat_res)
out_tif.SetGeoTransform(geotransform)
#定义投影
prj = osr.SpatialReference()
prj.ImportFromEPSG(4326)
out_tif.SetProjection(prj.ExportToWkt())
#数据导出
out_tif.GetRasterBand(1).WriteArray(tmp_arr[i]) #将数据写入内存,此时没有写入到硬盘
out_tif.FlushCache() #将数据写入到硬盘
out_tif = None #关闭tif文件
def main():
Input_folder = 'E:/桌面1/科创计划/cmip6/data'
Output_folder = 'E:/桌面1/科创计划/cmip6/nc2tif'
#读取所有数据
data_list = glob.glob(Input_folder+'*.nc')
for i in range(len(data_list)):
data = data_list[i]
nc2tif(data, Output_folder)
print(data+'转tif成功, 你真棒!')
main()