STMAS 数据可视化

import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.util import add_cyclic_point
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cbook as cbook
import matplotlib.image as image
from netCDF4 import Dataset
from cartopy.io.shapereader import Reader
import os
import numpy as np 
import pdb
import datetime,time
import pandas as pd
import xarray as xr
import math
import copy
from math import sin, asin, cos, radians, fabs, sqrt
import maskout
import geopandas as gpd
from matplotlib.font_manager import FontProperties
from matplotlib.transforms import offset_copy
import scipy
import scipy.ndimage as ndimage
from scipy.interpolate import griddata
from scipy.interpolate import Rbf
import xarray as xr
def interplot_data(data,lons,lats,glon,glat):

    data=(data.values).flatten()
    lats=(lats.values).flatten()
    lons=(lons.values).flatten()
    olon, olat = np.meshgrid(g_lon, g_lat)
    data = np.array(data).reshape(-1,1)
    lat = np.array(lats).reshape(-1,1)
    lon = np.array(lons).reshape(-1,1)
    points = np.concatenate([lon,lat],axis = 1)
    grid_data = griddata(points, data,(olon,olat), method='cubic')
    grid_data = grid_data[:,:,0]
    if olat[0,0]<olat[1,0]:
        olat = olat[-1::-1]
        grid_data = grid_data[-1::-1]
    return grid_data


if __name__ == '__main__':


    font = FontProperties(fname='/mnt/d/work/数据/SimHei.ttf', size=16)

    shp0='/mnt/d/work/数据/china_shp/cnhimap.shp'
    # shp1='/mnt/d/work/数据/浙江shp/Zj_County_Polygon1.shp'
    # shp1_geo=gpd.read_file(shp1)
    # shp1_geo.to_file('/mnt/d/work/数据/浙江shp/Zj_County_Polygon_utf-8.shp',encoding='utf-8')

    shp2='/mnt/d/work/数据/浙江shp/Zj_County_Polygon_utf-8.shp'
    shp222 = gpd.read_file(shp2)
    namelist='/mnt/d/work/20210106/data/stmas/namelist.txt'

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    mpl.rc('font', size=15, weight='normal')
    font = {'family': 'SimHei','weight': 'normal','size':40,}
    data_dir='/mnt/d/work/20210106/data/stmas/'
    for st_name in open(namelist):  
        st_name=st_name.strip('\n')
        print(st_name)  
        # st_iso=xr.open_dataset('./MSP3_PMSC_LAPS3KM_ME_L00_ZJ_202012291500_00000-00000.GR2', engine='cfgrib',
        # backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa'}})

        # st_sur=xr.open_dataset('./MSP3_PMSC_LAPS3KM_ME_L00_ZJ_202012291500_00000-00000.GR2', engine='cfgrib',
        # backend_kwargs={'filter_by_keys':{'typeOfLevel': 'surface'}})



        st_gr=xr.open_dataset(data_dir+st_name, engine='cfgrib',
        backend_kwargs={'filter_by_keys':{'typeOfLevel': 'heightAboveGround'}})

        st_iso=xr.open_dataset(data_dir+st_name, engine='cfgrib',
        backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa'}})

        time=st_gr.time.values+np.timedelta64(8, 'h')
        time=pd.to_datetime(str(time))
        timestr=str(time.year)+'年'+str(time.month)+'月'+str(time.day)+'时'+str(time.hour)+'时'+'（北京时）'

        u=st_iso['u'][2,:,:]
        v=st_iso['v'][2,:,:]
        temp=st_gr['t2m']-273.16
        lats=st_gr['latitude']
        lons=st_gr['longitude']

        nlat=110
        nlon=140
        g_lat=np.linspace(38,27,nlat)
        g_lon=np.linspace(110,124,nlon)

        u=interplot_data(u,lons,lats,g_lon,g_lat)
        v=interplot_data(v,lons,lats,g_lon,g_lat)
        temp=interplot_data(temp,lons,lats,g_lon,g_lat)
        # wsp=(u**2+v**2)**0.5
#        temp_s = scipy.interpolate.interp2d(g_lon,g_lat,temp,kind='linear')
        temp_s=ndimage.gaussian_filter(temp, sigma=3.0, order=0)

        olon, olat = np.meshgrid(g_lon, g_lat)


        proj= ccrs.PlateCarree() 
        fig = plt.figure(figsize=(10,8),dpi=300)  
        ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
        ax.add_geometries(Reader(shp0).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.3)
        ax.add_geometries(Reader(shp2).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='b', linewidth=0.3)
        #        ax.add_geometries(Reader(os.path.join(SHP_world, '10m_coastline.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.3)
        ax.set_xticks(np.arange(110,124,1))#需要显示的经度，一般可用np.arange
        ax.set_yticks(np.arange(27,40,1))#需要显示的纬度
        ax.xaxis.set_major_formatter(LongitudeFormatter())#将横坐标转换为经度格式
        ax.yaxis.set_major_formatter(LatitudeFormatter())#将纵坐标转换为纬度格式
        ax.tick_params(axis='both',labelsize=8,direction='in',length=2.75,width=0.55,right=True,top=True)#修改刻度样式
        ax.grid(linewidth=0.4, color='k', alpha=0.45, linestyle='--')#开启网格线
        extent=[118, 123, 27,31.5]
        ax.set_extent(extent,ccrs.PlateCarree())

        levels=np.arange(-10,11,1)
        tempf = ax.contourf(olon,olat,temp,levels=levels,cmap='coolwarm',transform=ccrs.PlateCarree())
        position=fig.add_axes([0.92,0.16,0.02,0.68])
        cb=fig.colorbar(tempf,cax=position,shrink=0.3,pad=0.01)

        # levels=np.arange(0,12,1)
        # wspf = ax.contourf(olon,olat,wsp_s,cmap='Blues',transform=ccrs.PlateCarree())
        levels=np.arange(0,1,1)
        temp_0 = ax.contour(olon,olat,temp_s,levels=levels,transform=ccrs.PlateCarree())

        barbs_inspace=2
        ax.barbs(olon[::barbs_inspace,::barbs_inspace],olat[::barbs_inspace,::barbs_inspace],
        u[::barbs_inspace,::barbs_inspace],v[::barbs_inspace,::barbs_inspace],
        linewidth=0.45,
        barb_increments={'half':2,'full':4,'flag':20},
        sizes=dict(spacing=0.1,width=0.05,emptybarb=0.02,height=0.35),length=3.5,zorder=5)

        ax.text(0.01,1.01,
                'STMAS 3KM Temperature(℃) & Wind',
                transform=ax.transAxes,
                fontsize=12)

        ax.text(0.67,1.01,
                timestr,
                transform=ax.transAxes,
                fontsize=14)

        logofile0='/mnt/d/work/logo/浙江服务中心-中.png'
        datafile0 = cbook.get_sample_data(logofile0,asfileobj=False)
        im = image.imread(datafile0)
        im[:, :, -1] = 1.0  # set the alpha channel
        fig.figimage(im, 360, 200, zorder=30)

        clip = maskout.shp2clip(tempf,
                                ax,
                                shpfile=shp2,
                                region='zhejiang',proj=proj)


        plt.savefig('./'+st_name+'.png',dpi=200)
        plt.close()