2022-11-09 10:16:06 +08:00
|
|
|
|
# 模块导入
|
|
|
|
|
|
import numpy as np
|
|
|
|
|
|
import netCDF4 as nc
|
|
|
|
|
|
from osgeo import gdal, osr, ogr
|
|
|
|
|
|
import os
|
2022-11-24 08:54:08 +08:00
|
|
|
|
# import glob
|
2022-11-09 10:16:06 +08:00
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
# 无论是日数据还是小时数居，变量名都是"SWR"
|
2022-11-10 16:17:29 +08:00
|
|
|
|
AOD_arr = np.asarray(nc_data_obj.variables["PAR"]) # 将SWR数据读取为数组
|
2022-11-09 10:16:06 +08:00
|
|
|
|
|
|
|
|
|
|
# 这个循环将所有Nodata的值（即-32768）全部改为0
|
|
|
|
|
|
for i in range(len(AOD_arr)):
|
|
|
|
|
|
for j in range(len(AOD_arr[0])):
|
|
|
|
|
|
if AOD_arr[i][j] == -32768:
|
|
|
|
|
|
AOD_arr[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)
|
|
|
|
|
|
|
|
|
|
|
|
# 此句代码必须有
|
|
|
|
|
|
AOD_arr = np.array([AOD_arr])
|
|
|
|
|
|
|
|
|
|
|
|
for i in range(len(AOD_arr[:])):
|
|
|
|
|
|
# 创建.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, 1, 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()) # 给新建图层赋予投影信息
|
|
|
|
|
|
|
|
|
|
|
|
# 数据写出
|
|
|
|
|
|
out_tif.GetRasterBand(1).WriteArray(AOD_arr[i]) # 将数据写入内存，此时没有写入硬盘
|
|
|
|
|
|
out_tif.FlushCache() # 将数据写入硬盘
|
|
|
|
|
|
out_tif = None # 注意必须关闭tif文件
|
