基于Python实现nc批量转tif格式

第基于Python实现nc批量转tif格式由于做项目需要运用到netCDF格式的气象数据，而ArcGIS中需要用栅格影像进行处理，对于较多的文件，ArcGIS一个个手动转换过于繁琐，因此我们采用Python进行转换，当然也可以采用matlab进行转换。

首先需要安装下面几个库：

importos

importnetCDF4asnc

importnumpyasnp

fromosgeoimportgdal,osr,ogr

importglob

我们可以在下面网址中寻找对应python安装版本的安装包，下载后，在pycharm控制台中直接安装即可。例如pipinstallnetCDF4-1.5.8-cp39-cp39-

win_amd64.whl

/~gohlke/pythonlibs/

安装之后即可进行转换：

defnc2tif(data,Output_folder):

tmp_data=nc.Dataset(data)#利用.Dataset()方法读取nc数据

print('tmp_data',tmp_data)

Lat_data=tmp_data.variables['lat'][:]

Lon_data=tmp_data.variables['lon'][:]

#print(Lat_data)

#print(Lon_data)

tmp_arr=np.asarray(tmp_data.variables['temp'])

#影像的左上角右下角坐标

Lonmin,Latmax,Lonmax,Latmin=[Lon_data.min(),Lat_data.max(),Lon_data.max(),Lat_data.min()]

#print(Lonmin,Latmax,Lonmax,Latmin)

#分辨率计算

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)

foriinrange(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)#WGS84

out_tif.SetProjection(prj.ExportToWkt())

#数据导出

out_tif.GetRasterBand(1).WriteArray(tmp_arr[i])#将数据写入内存，此时没有写入到硬盘

out_tif.FlushCache()#将数据写入到硬盘

out_tif=None#关闭tif文件

defmain():

Input_folder=r"E:\competition\航天宏图\2-meterairtemperature_CMFD\Data_forcing_01yr_010deg\\"

#Input_folder=r"E:\competition\航天宏图\2-meterairtemperature_CMFD\Data_forcing_01mo_010deg\\"

Output_folder=r"E:\competition\航天宏图\2-meterairtemperature_CMFD\tif\\"

#读取所有数据

data_list=glob.glob(os.path.join(Input_folder+'*.nc'))

print(data_list)

foriinrange(len(data_list)):

data=data_list[i]

nc2tif(data,Output_folder)

print(data+'转tif成功')

值得注意的是，tmp_arr=np.asarray(tmp_data.variables[temp])中的temp可以根据需要转换的波段进行选择，我们可以在读取数据之后print一下，找到对应的波段进行替换即可。如下图中我们应该选择的就是temp。

完成上述步骤即可得到所需的tif图像：

在上述代码中，经过处理的影像是倒置的，可能是处理过程中仿射矩阵读写错误导致的。因此我们可以在写入影像的时候，进行影像的垂直镜像操作即可：WriteArray(ndvi_arr_float[i][::-1])

defNC_to_tiffs(data,Output_folder):

nc_data_obj=nc.Dataset(data)

Lon=nc_data_obj.variables['lon'][:]

Lat=nc_data_obj.variables['lat'][:]

ndvi_arr=np.asarray(nc_data_obj.variables['temp'])

ndvi_arr_float=ndvi_arr.astype(float)/10000之间

#影像的左上角和右下角坐标

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)

foriinrange(len(ndvi_arr[:])):

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,N_Lon,N_Lat,1,gdal.GDT_Float32)

geotransform=(LonMin,Lon_Res,0,LatMax,0,-Lat_Res)

out_tif.SetGeoTransform(geotransform)

srs=osr.SpatialReference()

srs.ImportFromEPSG(4326)

out_tif.SetProjection(srs.ExportToWkt())

#数据写出

out_tif.GetRasterBand(1).WriteArray(ndvi_arr_float[i][::-1])#将数据写入内存，此时没有写入硬盘此处[::-1]用于图像的垂直镜像对称，避免图像颠倒

out_tif.FlushCache()#将数据写入硬盘

out_tif=None#注意必须关闭tif文件

这样便可以得到正确输出的影像：

当然，我们除了在写入时做垂直镜像操作之外，还可以利用python对图像进行几何变换来实现翻转。具体代码如下：

图像水平翻转：

#图像水平翻转

im_data_hor=np.flip(im_data,axis=2)

hor_path=train_image_path+"\\"+str(tran_num)+imageList[i][-4:]

writeTiff(im_data_hor,im_geotrans,im_proj,hor_path)

标签水平翻转：

#标签水平翻转

Hor=cv2.flip(label,1)

hor_path=train_label_path+"\\"+str(tran_num)+labelList[i][-4:]

cv2.imwrite(hor_path,Hor)

tran_num+=1

图像垂直翻转：

#图像垂直翻转

im_data_vec=np.flip(im_data,axis=1)

vec_path=train_image_path+"\\"+str(tran_num)+imageList[i][-4:]

writeTiff(im_data_vec,im_geotrans,im_proj,vec_path)

标签垂直翻转：

#标签垂直翻转

Vec=cv2.flip(label,0)

vec_path=train_label_path+"\\"+str(tran_num)+labelList[i][-4:]

cv2.imwrite(vec_path,Vec)

tran_num+=1

图像镜像翻转：

#图像对角镜像

im_data_dia=np.flip(im_data_vec,axis=2)

dia_path=train_image_path+"\\"+str(tran_num)+imageList[i][-4:]

writeTiff(im_data_dia,im_geotrans,im_proj,dia_pat