任务举例： 降水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()

任务举例： 如果转出的tif影像倒置。

操作： 加一行代码：

data = np.flipud(data)

任务举例： 若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)

编辑 | 南波婉帆 南波婉琳 南波婉琬

审核 | 南波婉琳

END 更多精彩关注GeoLab 219

​01准备工作：查看nc文件属性等。工具：Panoply、Matlab等软件。操作：1.使用Panoply 软件。2.使用Matlab软件。即可查看nc文件内各种属性；如果想单独查看变量：（以经度为例）02任务举例：降水nc数据转tif，坐标系为WGS 84。工具准备：Python GDAL和netCDF4环境。操作：替换代码中的属性变量名和路径即可。任务举例：降水nc数据转tif，坐标系为WGS 84。工具准备：Python GDAL和netC. 安装gdal库：pip install gdal 使用gdal库中的gdal_translate函数：gdal_translate -of G Tif f [input_file] [output_file] 如果需要，可以使用gdal库中的gdalwarp函数来重新定义输出 文件 的投影：gdalwarp -t_srs EPSG:43... from osgeo import gdal,osr var = 'SA' data = r'C:\Users\13290\Desktop\soil data\{}. nc '.format(var) f = nc .Dataset(data) var_lon = f['lon'][:] var_lat = f['lat'][:] data = f[var][0, :] data_arr = np.asarray(data) x = nc _file.variables['x'][:] y = nc _file.variables['y'][:] transform = (x[0], x[1]-x[0], 0, y[0], 0, y[1]-y[0]) # 创建 tif 文件 driver = gdal.GetDriverByName('G Tif f') tif _file = driver.Create('your_ tif _file. tif ', data.shape[1], data.shape[0], 1, gdal.GDT_Float32) # 设置 tif 文件 的元数据 tif _file.SetGeoTransform(transform) tif _file.SetProjection( nc _file.variables['your_variable'].get nc attr('crs')) # 将数据写入 tif 文件 tif _file.GetRasterBand(1).WriteArray(data) # 关闭 文件 nc _file.close() tif _file = None 注意，以上 代码 仅供参考，具体操作需要根据实际情况进行调整。