制作手机网站哪家好多用户商城源码下载
本文介绍基于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