??本文介紹基于Python中GDAL模塊，批量將大量.nc格式的柵格文件轉換為.tif格式，并解決可能出現的轉換后圖像顛倒、鏡像、翻轉等問題。

??最近，需要批量將大量.nc格式的柵格文件轉換為.tif格式。如下圖所示，有多個待轉換的.nc格式文件，且對于每一個.nc格式文件，其都含有多個時相的數據。

??其實，對于.nc格式轉.tif格式，除了gdal庫之外，還有其他非常多成熟、方便的Python庫或R語言庫可以實現，且這些庫都比gdal庫用起來方便——甚至安裝過程也比gdal庫方便；但是，我在一開始用這些其他庫進行格式轉換時發現：對于我的.nc格式數據，若用其他庫進行轉換，最終得到的.tif圖像并不是正確的——其要么是行列數（也就是經緯度）都反了，要么是經度或緯度其中一個是反著的；包括Python的netCDF4庫、xarray庫、rioxarray庫，以及R語言的ncdf庫等，都存在這個問題。

??如下圖所示，這個就是我遇到的其中一種情況。可以看到，這個結果數據的經度倒是沒錯，但是緯度搞反了，所以全球的圖像是反著來的，南極跑到北極去了。

??針對這種情況，我大致看了一下原本的.nc格式數據，感覺問題應該就是出在.nc數據上——比如可能對于一些學者自行生產的科學數據產品，其在數據生產完畢、封裝為.nc格式時，經緯度是反著寫的；而上述提及的那些成熟的.nc格式轉.tif格式的庫，他們默認是正著解析經度和緯度的，所以出現了上述問題。

??所以，為了解決這個問題，那就不太好用這些封裝好的.nc格式轉.tif格式庫了，而是需要將.nc格式數據直接提取為類似于array格式的矩陣數據，然后手動進行矩陣變換，再將其導出為.tif——那這樣的話，就只有gdal庫符合要求了。

??本文所用代碼如下。

import os
import numpy as np
import netCDF4 as nc
from osgeo import gdal
from osgeo import osrnc_folder = R"e:\LTDR\N07"
tif_folder = R"d:\99_Temp\NDVI\SPEI\TIFF"
data_name = "spei"
nodata = -9999
scale = 10000
res = 0.25if not os.path.exists(tif_folder):os.makedirs(tif_folder)driver = gdal.GetDriverByName('GTiff')for nc_file in os.listdir(nc_folder):if nc_file.endswith(".nc"):nc_file_path = os.path.join(nc_folder, nc_file)year = nc_file.split("_")[2]with nc.Dataset(nc_file_path) as file:data = np.array(file[data_name][:])times = np.array(file['time'][:])lats = np.array(file['lat'][:])lons = np.array(file['lon'][:])lat_min, lat_max = lats.min(), lats.max()lon_min, lon_max = lons.min(), lons.max()lat_num = len(lats)lon_num = len(lons)lat_res = reslon_res = reslat_origin = lat_max + lat_res / 2lon_origin = lon_min - lon_res / 2for time in range(len(times)):daily_data = data[time, :, :]# scale the data if necessary# daily_data = np.round(data[time, :, :] * scale).astype(np.int32)# Flip the data to match the expected orientation, here are three options:daily_data = np.flipud(daily_data.T)daily_data = np.flipud(daily_data)daily_data = daily_data.Toutput_tif_path = os.path.join(tif_folder,f"SPEI_{year}{str(time + 1).zfill(3)}.tif")outRaster = driver.Create(output_tif_path, lon_num, lat_num, 1, gdal.GDT_Float32)# outRaster = driver.Create(output_tif_path, lon_num, lat_num, 1, gdal.GDT_Int32)outRaster.SetGeoTransform([lon_origin, lon_res, 0, lat_origin, 0, -lat_res])sr = osr.SpatialReference()sr.SetWellKnownGeogCS('WGS84')outRaster.SetProjection(sr.ExportToWkt())band = outRaster.GetRasterBand(1)band.SetNoDataValue(nodata)# scale the data if necessary# band.SetNoDataValue(nodata * scale)band.WriteArray(daily_data)print(output_tif_path, ' finished')del outRaster

??其中，nc_folder就是.nc格式文件所在目錄，tif_folder是輸出.tif格式文件的目錄，data_name是要提取的變量名，nodata表示填充值，scale是縮放系數（例如假設原本的.nc格式文件是浮點數，乘以10000并取整后，存為整數可節省空間），res則是空間分辨率。

??代碼整體思路也很簡單，主要就是需要關注daily_data = np.flipud(daily_data.T)這句代碼以及其下方的2行代碼。這里就是將原本.nc格式文件數據加以變換的地方，這里列出了3種變換方法，分別為先轉置、后上下顛倒，以及直接上下顛倒，還有直接轉置。這里之所以列出3種方法，是因為我當時要轉換的.nc格式數據產品有很多種，為了方便就將不同種對我有效果的變換方法都寫上了；大家使用代碼時，需要注意選擇自己適合的變換方法。還有一種情況，就是可能圖像還會出現左右顛倒的問題，也就是緯度沒問題、但經度反了——不過這種情況感覺一般不會遇到，所以當時就沒寫變換的代碼；如果大家遇到了，那就需要額外對array進行左右變換。

??至此，大功告成。

歡迎關注：瘋狂學習GIS