01
准备工作：查看nc文件属性等。
工具：Panoply、Matlab等软件。
2.使用Matlab软件。
即可查看nc文件内各种属性；
如果想单独查看变量：（以经度为例）
02
任务举例：降水nc数据转tif，坐标系为WGS 84。
工具准备：Python GDAL和netCDF4环境。
操作：替换代码中的属性变量名和路径即可。
任务举例：降水nc数据转tif，坐标系为WGS 84。
工具准备：Python GDAL和netCDF4环境。
操作：替换代码中的属性变量名和路径即可。
-*- coding: UTF-8 -*-
import os
import netCDF4 as nc
import numpy as np
from osgeo import gdal, osr, ogr
import glob
def nc2tif(data, Output_folder):
pre_data = nc.Dataset(data) # 利用.Dataset()读取nc数据
Lat_data = pre_data.variables['lat'][:]
# print(Lat_data)
Lon_data = pre_data.variables['lon'][:]
# print(Lon_data)
pre_arr= np.asarray(pre_data.variables['pre'])#属性变量名
# 影像的左上角&右下角坐标
Lonmin, Latmax, Lonmax, Latmin = [Lon_data.min(), Lat_data.max(), Lon_data.max(), Lat_data.min()]
# Lonmin, Latmax, Lonmax, Latmin
# 分辨率计算
Num_lat = len(Lat_data)
Num_lon = len(Lon_data)
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)
for i in range(len(pre_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(pre_arr[i]) # 将数据写入内存
out_tif.FlushCache() # 将数据写入到硬盘
out_tif = None # 关闭tif文件
def main():
Input_folder = r'G:/pre_2020/'
Output_folder = r'G:/pre_2020/'
data_list = glob.glob(os.path.join(Input_folder, '*.nc'))
data_list.sort()
# 读取所有数据
# data_list = glob.glob(Input_folder + '*.nc')
print(data_list)
for i in range(len(data_list)):
data = data_list[i]
nc2tif(data, Output_folder)
print(data + '转tif成功')
main()
03
任务举例：如果转出的tif影像倒置。
操作：加一行代码：
data = np.flipud(data)
示例：
04
任务举例：若nc文件范围跨了180度经线，无法确定具体范围和分辨率；或nc文件的坐标系并不是WGS84。
操作：利用Arcgis转出一幅范围和投影正确的栅格图像，再用代码。（Arcgis无法将一个nc文件导出为不同时间的tif文件）。
步骤一：利用arcgis将nc转tif。
系统工具箱-Multidimension Tools-创建NetCDF栅格图层。
成功后右击图层-数据-导出数据，即可导出栅格数据集。
步骤二：修改路径和变量名，在python中运行下列代码即可。
-*- coding: UTF-8 -*-
import os
import netCDF4 as nc
import numpy as np
from osgeo import gdal, osr, ogr
import glob
def WriteTiff(im_data, inputdir, path):
raster = gdal.Open(inputdir)
im_width = raster.RasterXSize # 栅格矩阵的列数
im_height = raster.RasterYSize # 栅格矩阵的行数
im_bands = raster.RasterCount # 波段数
im_geotrans = raster.GetGeoTransform() # 获取仿射矩阵信息
im_proj = raster.GetProjection() # 获取投影信息
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])
im_bands, im_height, im_width = im_data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, im_width, im_height, im_bands, datatype)
#print(dataset)
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
infile = "D:/tmax/daymet_v4_daily_na_tmax_1980.nc" #nc文件路径
data_set = nc.Dataset(infile) # 读取nc文件信息
time = data_set.variables["time_bnds"][:] #获取时间序列
读取样本tif文件的地理信息
intif = "D:/tmax_Layer1.tif" #标准的tif文件路径
for i in range(0, len(time)):
value_data = data_set.variables['tmax'][i]
# 将缺失值改为0
data = value_data.data
mask = value_data.mask
data[np.where(mask == True)] = 0
outputname = "D:/tmax2/" + 'tmax' + str(1980)+str(i+1).zfill(3) + ".tif" #输出文件的命名方式
WriteTiff(data, intif, outputname)
print(outputname)
undefined
编辑 | 南波婉帆 南波婉琳 南波婉琬
审核 | 南波婉琳
END 更多精彩关注GeoLab 219
Original: https://blog.csdn.net/weixin_42776126/article/details/124177572
Author: GeoLab 219
Title: 代码分享 Python | 批量nc文件转tif
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/632543/
转载文章受原作者版权保护。转载请注明原作者出处！