读取智能网格数据并叠加交通线路和格点值

专供浙江省交警总队最低气温.png

代码如下：

#读ec数据并画图hgt,wind,rh
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 matplotlib.transforms import offset_copy
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 maskout
import geopandas as gpd
if __name__ == '__main__':
    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,}


    file_dir='/mnt/d/work/other/liangxn/7/data/'
    file_name='Ruler.20201231080000.024.TMin24.nc'
    init_time='起报时间：2020-12-29 20:00' 
    title_txt='2021年1月2日24小时最低气温'
    client_name='专供浙江省交警总队'

    logofile='/mnt/d/work/浙江省气象服务中心.png'
    shp0='/mnt/d/work/数据/浙江Shp/Zj_County_Polygon_utf-8.shp'
    shp1='/mnt/d/work/数据/浙江道路/0.浙江/杭州市道路数据/高速公路.shp'
    shp2='/mnt/d/work/数据/浙江道路/0.浙江/湖州市道路数据/高速公路.shp'
    shp3='/mnt/d/work/数据/浙江道路/0.浙江/嘉兴市道路数据/高速公路.shp'
    shp4='/mnt/d/work/数据/浙江道路/0.浙江/金华市道路数据/高速公路.shp'
    shp5='/mnt/d/work/数据/浙江道路/0.浙江/丽水市道路数据/高速公路.shp'
    shp6='/mnt/d/work/数据/浙江道路/0.浙江/宁波市道路数据/高速公路.shp'
    shp7='/mnt/d/work/数据/浙江道路/0.浙江/衢州市道路数据/高速公路.shp'
    shp8='/mnt/d/work/数据/浙江道路/0.浙江/绍兴市道路数据/高速公路.shp'
    shp9='/mnt/d/work/数据/浙江道路/0.浙江/台州市道路数据/高速公路.shp'
    shp10='/mnt/d/work/数据/浙江道路/0.浙江/温州市道路数据/高速公路.shp'
    shp11='/mnt/d/work/数据/浙江道路/0.浙江/舟山市道路数据/高速公路.shp'

    proj= ccrs.PlateCarree() 
    fig = plt.figure(figsize=(10,8),dpi=300)  
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
    traffic_edgecolor='r'
    traffic_linewidth=0.5
    ax.add_geometries(Reader(shp0).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.1)
    ax.add_geometries(Reader(shp1).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp2).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp3).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp4).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp5).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp6).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp7).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp8).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp9).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp10).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)
    ax.add_geometries(Reader(shp11).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor=traffic_edgecolor, linewidth=traffic_linewidth)


    nc = xr.open_dataset(file_dir+file_name)
    #nc = xr.open_dataset(u'/mnt/d/work/other/liangxn/7/20210102/privonce/Ruler.20210102200000.024.TMin24.nc')
    lons = nc['lon']
    lats = nc['lat']
    tmin=nc['TMin24']
    olon,olat=np.meshgrid(lons,lats)


    ax.set_xticks([118,118.5,119,119.5,120,120.5,121,121.5,122,122.5,123])#需要显示的经度，一般可用np.arange
    ax.set_yticks([27.5,28,28.5,29,29.5,30,30.5,31,31.5,32])#需要显示的纬度
    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, 26.9,31.5]
    ax.set_extent(extent,ccrs.PlateCarree())

    print(tmin.max(),tmin.min())
    levels=np.arange(-18,0,2)
    cf = ax.contourf(olon,olat,tmin,levels=levels,cmap='Blues_r',transform=ccrs.PlateCarree())
    position=fig.add_axes([0.92,0.16,0.02,0.68])
    cb=fig.colorbar(cf,cax=position,shrink=0.3,pad=0.01)
#    ax.set_title('起报时间：2020-12-26 20:00        单位：℃         专供浙江省高速交警总队',fontsize=12)
    geodetic_transform = ccrs.Geodetic()._as_mpl_transform(ax)
    text_transform = offset_copy(geodetic_transform, units='dots', x=-25)
    inter=5

    for i in range(25,len(lons)-25,inter):
        for j in range(20,len(lats)-10,inter):
            ax.text(lons[i],
                    lats[j],
                    int(tmin[j,i].values),
                    verticalalignment='top',
                    transform=text_transform,
                    fontsize=10)

    ax.text(0.01,1.01,
            init_time,
            transform=ax.transAxes,
            fontsize=12)

    ax.text(0.86,1.01,
            '单位：℃',
            transform=ax.transAxes,
            fontsize=12)

    ax.text(0.65,0.95,
            client_name,
            transform=ax.transAxes,
            fontsize=12)
    # 白化

#     zj = gpd.read_file(shp0) 
#     zj.to_file('/mnt/d/work/数据/浙江边界省市县/行政境界-市-面.shp',encoding='utf8')  

    # 白化
    clip = maskout.shp2clip(cf,
                            ax,
                            shpfile=shp0,
                            region='zhejiang',proj=proj)

    datafile = cbook.get_sample_data(logofile,asfileobj=False)
    im = image.imread(datafile)
    im[:, :, -1] = 1.0  # set the alpha channel
    fig.figimage(im, 550, 280, zorder=22)
    fig.suptitle(title_txt)
    plt.savefig('./'+client_name+'最低气温.png')
    plt.close()