import numpy as np
import netCDF4 as nc
from osgeo import gdal, osr, ogr
import os
import glob
def NC_to_tiffs(data, Output_folder):
# 读取一下基本信息
nc_data_obj = nc.Dataset(data)
Lon = nc_data_obj.variables["longitude"][:]
Lat = nc_data_obj.variables["latitude"][:]
# 读取变量的时候,会自动根据scale factor对数值进行还原,但是Nodata的栅格会存储为-32768
# 无论是日数据还是小时数居,变量名都是"每个波段的数据"
band1 = np.asarray(nc_data_obj.variables["albedo_01"]) # 将数据读取为数组
band2 = np.asarray(nc_data_obj.variables["albedo_02"])
band3 = np.asarray(nc_data_obj.variables["albedo_03"])
band4 = np.asarray(nc_data_obj.variables["albedo_04"])
# 这个循环将所有Nodata的值(即-32768)全部改为0
for i in range(len(band1)):
for j in range(len(band1[0])):
if band1[i][j] == -32768:
band1[i][j] = 0.0
# 影像的四秩
LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]
# 分辨率计算,其实可以写死,都是2401*2401
N_Lat = len(Lat)
N_Lon = len(Lon)
Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
bands = [band1, band2, band3, band4]
# 创建.tif文件
driver = gdal.GetDriverByName("GTiff")
TIFF_name = os.path.basename(data)
out_tif_name = Output_folder + "\\" + TIFF_name.split("_")[1] + "_" + TIFF_name.split("_")[2] + ".tif"
out_tif = driver.Create(out_tif_name, N_Lon, N_Lat, 4, gdal.GDT_Float32)
# 设置影像的显示范围
# -Lat_Res一定要是负的
geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
out_tif.SetGeoTransform(geotransform)
# 获取地理坐标系统信息,用于选取需要的地理坐标系统
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84",AUTHORITY["EPSG","4326"]
out_tif.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息
for index, band in enumerate(bands):
# 此句代码必须有
band = np.array([band])
for i in range(len(band[:])):
# 数据写出
out_tif.GetRasterBand(index + 1).WriteArray(band[i]) # 将每个波段的数据写入内存,此时没有写入硬盘
out_tif.FlushCache() # 最终将数据写入硬盘
out_tif = None # 注意必须关闭tif文件
if __name__ == "__main__":
# Download_Path用于存储下载的原始数据
Download_Path = "D:/data/Himawari-8"
# Analysis_Path用于存储处理后的数据(即转为TIFF后的数据)的文件夹
Analysis_Path = "D:/data/Himawari-8/tiff_file"
# 下面开始数据处理
# 读取所有nc数据
data_list = glob.glob(Download_Path + "\\*.nc")
# for循环完成解析
for i in range(len(data_list)):
data = data_list[i]
NC_to_tiffs(data, Analysis_Path)
print(data + "-----转tif成功")
print("----转换结束----")
葵花8号全球影像nc产品按波段转成tiff
最新推荐文章于 2023-08-16 16:21:16 发布