# 导入包
from netCDF4 import Dataset
import numpy as np
from osgeo import gdal,osr
import os
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____20221029T060755_20221029T061255_26131_03_020401_20221029T065050.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.keys())
# print (dataset_nc.groups['PRODUCT'].variables['sulfurdioxide_total_vertical_column_precision'])
#获取so2对象
lons =dataset_nc.groups['PRODUCT'].variables['longitude'][:][0,:,:]
lats = dataset_nc.groups['PRODUCT'].variables['latitude'][:][0,:,:]
no2 = dataset_nc.groups['PRODUCT'].variables['sulfurdioxide_total_vertical_column_precision'][0,:,:]
print (lons.shape)
print (lats.shape)
print (no2.shape)


no2_units = dataset_nc.groups['PRODUCT'].variables['sulfurdioxide_total_vertical_column_precision'].units
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from mpl_toolkits.basemap import Basemap
lon_0 = lons.mean()
lat_0 = lats.mean()

m = Basemap(width=5000000,height=3500000,
            resolution='l',projection='stere',\
            lat_ts=40,lat_0=lat_0,lon_0=lon_0)

xi, yi = m(lons, lats)

# Plot Data
cs = m.pcolor(xi,yi,np.squeeze(no2),norm=LogNorm(), cmap='jet')

# Add Grid Lines
m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries
m.drawcoastlines()
m.drawstates()
m.drawcountries()

# Add Colorbar
cbar = m.colorbar(cs, location='bottom', pad="10%")
cbar.set_label(no2_units)

# Add Title
plt.title('SO2 in atmosphere(2022.10.29)')
plt.show()