# _*_ coding:utf-8 _*_
#----------------------导入数据读取和处理的模块-----------------------
from netCDF4 import Dataset
from pathlib import Path
import glob
import xarray as xr
import numpy as np
import pandas as pd
from wrf import getvar, ALL_TIMES, latlon_coords, xy_to_ll, ll_to_xy, smooth2d, \
get_cartopy,cartopy_xlim,cartopy_ylim,to_np,to_xy_coords,interplevel,CoordPair,vertcross
#----------------------导入画图相关函数-------------------------------
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import matplotlib.ticker as ticker
from cartopy import mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import cartopy.io.shapereader as shapereader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter,LongitudeLocator,LatitudeLocator
from pylab import *
#-------------------- 导入颜色包---------------------------------------
import seaborn as sns
from matplotlib import cm
from matplotlib.colors import ListedColormap
#------------------- 导入插值模块----------------------------------------
from scipy.interpolate import interp1d #引入scipy中的一维插值库
from scipy.interpolate import griddata #引入scipy中的二维插值库
from scipy.interpolate import interp2d
#----------------------字体设置--------------------------------------------
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 15
#输入路径,读取nc数据属性
path = "J:/0hsyPaper202406/4data/4-wrfout-231106/copy_merge/wrfout_wuwei1.nc"
ds=xr.open_dataset(path)
ds.dims
#nc数据coords变量名重命名。
df = ds.rename({"Time":"time"})
df.dims
dft2=df["T2"]
dft2["XTIME"].data
# dft2["XLAT"].data
t2=dft2.loc[29:,9:108,9:88]
t2["XTIME"].data
t2_spr = t2.loc[0:92]
# t2_spr["XTIME"].data
t2_sum = t2.loc[92:184]
# t2_sum["XTIME"].data
t2_aut = t2.loc[184:275]
# t2_win["XTIME"].data
t2_spr_meant = t2_spr.mean(dim="time")-273.15 #(-14.80969238, 9.49023438)
print(t2_spr_meant)
t2_sum_meant = t2_sum.mean(dim="time")-273.15 #(-3.60955811, 22.57196045)
t2_aut_meant = t2_aut.mean(dim="time")-273.15 #(-14.19775391, 6.83230591)
lats,lons=latlon_coords(t2)
# lats,lons
fig=plt.figure(figsize=(20,18))
cart_proj = ccrs.LambertConformal(central_longitude=103, standard_parallels=(30, 60))
###################################################################################################
#t2_spr_meant
###################################################################################################
t2min=-15
t2max=25
ax1=plt.subplot(1,3,1,projection=cart_proj)
#1、图题########################
ax1.set_title("Spring T2 temperature",fontsize=20)
ax1.set_title('unit: ℃',loc='right',fontsize=16)
#2、叠加区域shp#######################################
shp="F:/市.shp"
cnmap=cfeat.ShapelyFeature(shapereader.Reader(shp).geometries(), ccrs.PlateCarree(), edgecolor="black", facecolor='none')
ax1.add_feature(cnmap, linewidth=1)
#3、设置模拟区域范围与网格形式########################
levels1=np.arange(-15,15)
im1=plt.pcolormesh(to_np(lons[0]), to_np(lats[0]), to_np(t2_spr_meant), vmin=t2min, vmax=t2max,
cmap=get_cmap("coolwarm"),transform=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax1.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lon_format=LongitudeFormatter()
lat_format=LatitudeFormatter()
gl1 = ax1.gridlines(draw_labels = True, linestyle=":", crs=ccrs.PlateCarree(), linewidth=1, color="grey",
x_inline=False, y_inline=False)
#gl.top_labels = False #关闭上部经纬度标签
gl1.right_labels = False
gl1.top_labels = False
# gl.xlocator = mticker.FixedLocator(np.arange(70,140,10))
# gl.ylocator = mticker.FixedLocator(np.arange(10,60,10))
gl1.xlabel_style ={"size":18, "color":"black"} #修改经纬度坐标网格字体大小
gl1.ylabel_style ={"size":18, "color":"black"} #修改经纬度坐标网格字体大小
gl1.rotate_labels = False
#4、设置colorbar###########################
# cbar1=plt.colorbar(ax=ax1,shrink=0.98)
###################################################################################################
#t2_sum_meant
###################################################################################################
ax2=plt.subplot(1,3,2,projection=cart_proj)
#1、图题########################
ax2.set_title("Summer T2 temperature",fontsize=20)
ax2.set_title('unit: ℃',loc='right',fontsize=16)
#2、叠加区域shp#######################################
shp="F:/市.shp"
cnmap=cfeat.ShapelyFeature(shapereader.Reader(shp).geometries(), ccrs.PlateCarree(), edgecolor="black", facecolor='none')
ax2.add_feature(cnmap, linewidth=1)
#3、设置模拟区域范围与网格形式########################
levels2=np.arange(0,25)
im2=plt.pcolormesh(to_np(lons[0]), to_np(lats[0]), to_np(t2_sum_meant), vmin=t2min, vmax=t2max,
cmap=get_cmap("coolwarm"),transform=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax2.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lon_format=LongitudeFormatter()
lat_format=LatitudeFormatter()
gl2 = ax2.gridlines(draw_labels = True, linestyle=":", crs=ccrs.PlateCarree(), linewidth=1, color="grey",
x_inline=False, y_inline=False)
#gl.top_labels = False #关闭上部经纬度标签
gl2.right_labels = False
gl2.top_labels = False
# gl.xlocator = mticker.FixedLocator(np.arange(70,140,10))
# gl.ylocator = mticker.FixedLocator(np.arange(10,60,10))
gl2.xlabel_style ={"size":18, "color":"black"} #修改经纬度坐标网格字体大小
gl2.ylabel_style ={"size":18, "color":"black"} #修改经纬度坐标网格字体大小
gl2.rotate_labels = False
#4、设置colorbar###########################
# cbar2=plt.colorbar(ax=ax2,shrink=0.98)
###################################################################################################
#t2_win_meant
###################################################################################################
ax3=plt.subplot(1,3,3,projection=cart_proj)
#1、图题########################
ax3.set_title("Autumn T2 temperature",fontsize=20)
ax3.set_title('unit: ℃',loc='right',fontsize=16)
#2、叠加区域shp#######################################
shp="F:/市.shp"
cnmap=cfeat.ShapelyFeature(shapereader.Reader(shp).geometries(), ccrs.PlateCarree(), edgecolor="black", facecolor='none')
ax3.add_feature(cnmap, linewidth=1)
#3、设置模拟区域范围与网格形式########################
levels3=np.arange(-15,25)
im3=plt.pcolormesh(to_np(lons[0]), to_np(lats[0]), to_np(t2_aut_meant), vmin=t2min, vmax=t2max,
cmap=get_cmap("coolwarm"),transform=ccrs.PlateCarree())
ax3.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax3.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lon_format=LongitudeFormatter()
lat_format=LatitudeFormatter()
gl3 = ax3.gridlines(draw_labels = True, linestyle=":", crs=ccrs.PlateCarree(), linewidth=1, color="grey",
x_inline=False, y_inline=False)
#gl.top_labels = False #关闭上部经纬度标签
gl3.right_labels = False
gl3.top_labels = False
# gl.xlocator = mticker.FixedLocator(np.arange(70,140,10))
# gl.ylocator = mticker.FixedLocator(np.arange(10,60,10))
gl3.xlabel_style ={"size":18, "color":"black"} #修改经纬度坐标网格字体大小
gl3.ylabel_style ={"size":18, "color":"black"} #修改经纬度坐标网格字体大小
gl3.rotate_labels = False
#4、设置colorbar###########################
# cbar3=plt.colorbar(ax=ax3,shrink=0.98)
#####################################################################################################################
fig.subplots_adjust(right=0.9) #前面三个子图的总宽度为全部宽度的0.9,剩下的0.1用来放置colorbar
#colorbar左下宽高
l=0.92
b=0.31
w=0.018
h=1-2*b
#对应l,b,w,h,设置colorbar位置
rect = [l,b,w,h]
cbar_ax = fig.add_axes(rect)
cb=plt.colorbar(im3, cax=cbar_ax, shrink=0.68)
#设置colorbar标签字体等
cb.ax.tick_params(labelsize=16)
font = {"family":"serif",
"color": "black",
"weight":"normal"}
# cb.set_label("T", fontdict=font, loc="top")
cb.ax.set_title("T2 ℃")
plt.savefig("t2_season_mean_wuwei1_pcolormesh.png", dpi=900,bbox_inches = 'tight')