# 导入包
from netCDF4 import Dataset
import numpy as np
from osgeo import gdal,osr
import os
import math
os.environ['CPL_ZIP_ENCODING'] = 'UTF-8'
os.environ['PROJ_LIB'] = r'D:\anaconda3\envs\pytorch_GPU\Lib\site-packages\pyproj\proj_dir\share\proj'
os.environ['GDAL_DATA'] = r'D:\anaconda3\envs\pytorch_GPU\Library\share'

# nc文件路径
Nc_file_path  = r"C:\Users\1\Desktop\S5P_NRTI_L2__SO2____20221102T063255_20221102T063755_26188_03_020401_20221102T071551.nc"
# 读取nc文件
dataset_nc  = Dataset(Nc_file_path , mode='r')
# 打印 nc以查看文件是否已加载及nc文件存储结构及信息
# print (dataset_nc)
# 打印dataset_nc的groups属性
# print (dataset_nc.groups)
# 打印dataset_nc的groups属性中的PRODUCT
print(dataset_nc.groups['PRODUCT'].variables['corner'])
nc_so2= dataset_nc.groups['PRODUCT'].variables['sulfurdioxide_total_vertical_column_precision']
#获取so2对象
so2_value = dataset_nc.groups['PRODUCT'].variables['sulfurdioxide_total_vertical_column_precision'][0,:,:]
#无效值
_FillValue = nc_so2._FillValue
#数据的单位
so2_pixel = nc_so2.units
so2_arr= np.asarray(so2_value).astype(float)

# a=so2_arr.astype(np.uint8)
so2_arr =np.where(so2_arr!=9.96921e+36,so2_arr,0)
so2_arry_max = so2_arr.max()
# so2_arr2 =np.where(so2_arr1!=-9999,so2_arr1 *10,0)
# 获取其数据对应的经纬度,Sentinel-5P 数据输出是一个 3D 数组，因此我们需要使用 [0，：，：] 代码将其切片为 2D 数组
# 经度
Lon_data = dataset_nc.groups['PRODUCT'].variables['longitude'][0,:,:]
# 纬度
Lat_data = dataset_nc.groups['PRODUCT'].variables['latitude'][0,:,:]
# -------------------进行输出tif文件的准备---------------------
# 影像的左上角&右下角坐标
Lonmin, Latmax, Lonmax, Latmin = [Lon_data.min(), Lat_data.max(),Lon_data.max(),Lat_data.min()]
# Lonmin, Latmax, Lonmax, Latmin
# 分辨率计
# w
Num_lat = so2_arr.shape[0]
# h
Num_lon = so2_arr.shape[1]
Lat_res = (Latmax - Latmin) / (float(Num_lat) - 1.)
Lon_res = (Lonmax - Lonmin) / (float(Num_lon) - 1.)



# 创建tif文件
driver = gdal.GetDriverByName('GTiff')
# 输出名称
out_tif_name = r'C:\Users\1\Desktop\testso11.tif'
#使用gdal驱动器，创建tiff文件
out_tif = driver.Create(out_tif_name,  Num_lat, Num_lon ,1, gdal.GDT_Float32)

# 设置影像的显示范围
# 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_shape = out_tif.GetRasterBand(1).shape
out_tif.GetRasterBand(1).WriteArray(so2_arr.T)  # 将数据写入内存
out_tif.FlushCache()  # 将数据写入到硬盘
out_tif = None  # 关闭tif文件