- 积分
- 970
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-10-27
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2020-12-30 22:12:29
|
显示全部楼层
# -*- coding: utf-8 -*-
"""
读取NC文件,变量有位势高度,U风,V风
"""
import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
from matplotlib import pyplot as plt
from netCDF4 import Dataset
import matplotlib as mpl
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
mpl.rcParams['font.sans-serif'] = ['SimHei'] #设置简黑字体
mpl.rcParams['axes.unicode_minus'] = False #正常显示负号
#-----读取数据
nc_obj2014=Dataset('./data/2014.nc')
nc_obj2015=Dataset('./data/2015.nc')
nc_obj2016=Dataset('./data/2016.nc')
nc_obj2017=Dataset('./data/2017.nc')
nc_obj2018=Dataset('./data/2018.nc')
nc_obj2019=Dataset('./data/2019.nc')
lat = (nc_obj2014.variables['latitude'][:])
lon = (nc_obj2014.variables['longitude'][:])
level = (nc_obj2014.variables['level'][:])
lons,lats = np.meshgrid(lon,lat)
time2015 = (nc_obj2015.variables['time'][:])
hgt2014 = (nc_obj2014.variables['z'][:,9,:,:])*0.1
hgt2015 = (nc_obj2015.variables['z'][:,9,:,:])*0.1
hgt2016 = (nc_obj2016.variables['z'][:,9,:,:])*0.1
hgt2017 = (nc_obj2017.variables['z'][:,9,:,:])*0.1
hgt2018 = (nc_obj2018.variables['z'][:,9,:,:])*0.1
hgt2019 = (nc_obj2019.variables['z'][:,9,:,:])*0.1
#----绘图
fontNum=12
contourline = range(11970,12600,30)
region = [100,120,15,25]
proj=ccrs.PlateCarree()
fig = plt.figure(figsize=(6,6),dpi=600)
#----第一张图
ax1 = fig.add_axes([0.06,0.62,0.45,0.4],
projection=ccrs.LambertConformal(central_longitude=110.0,
central_latitude=20,standard_parallels=(15,25)))
ax1.set_extent(region,crs=ccrs.PlateCarree())
ax1.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'],alpha=0.5)
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.STATES.with_scale('50m'))
cs1 = ax1.contour(lons[280:361,160:321],lats[280:361,160:321],hgt2015[10,280:361,160:321],
transform=proj,levels=contourline,linewidths=1,colors="k")
cs1.clabel(contourline[0::2],fontsize=fontNum,colors="k",fmt='%.f')
gl = ax1.gridlines(crs=proj,ylocs=np.arange(region[2],region[3]+5,5),
xlocs=np.arange(region[0],region[1]+5,5),draw_labels=True)
ax1.yaxis.set_ticks_position("right")
#ax1.set_xticks(range(100,120,5),crs=proj)
#ax1.set_yticks(range(15,25,5),crs=proj)
ax1.set_title("201511",fontsize=15)
#----第二张图
ax2 = fig.add_axes([0.52,0.62,0.45,0.4],
projection=ccrs.LambertConformal(central_longitude=110.0,
central_latitude=20,standard_parallels=(15,25)))
ax2.set_extent(region,crs=ccrs.PlateCarree())
cs2 = ax2.contour(lons[280:361,160:321],lats[280:361,160:321],hgt2015[11,280:361,160:321],
transform=proj,levels=contourline,linewidths=1,colors="k")
cs2.clabel(contourline[0::2],fontsize=fontNum,colors="k",fmt='%.f')
ax2.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'],alpha=0.5)
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax2.add_feature(cfeature.STATES.with_scale('50m'))
gl2 = ax2.gridlines(crs=proj,ylocs=np.arange(region[2],region[3]+5,5),
xlocs=np.arange(region[0],region[1]+5,5),draw_labels=False)
ax2.set_title("201512",fontsize=15)
#----第三张图
ax3 = fig.add_axes([0.06,0.32,0.45,0.4],
projection=ccrs.LambertConformal(central_longitude=110.0,
central_latitude=20,standard_parallels=(15,25)))
ax3.set_extent(region,crs=ccrs.PlateCarree())
cs3 = ax3.contour(lons[280:361,160:321],lats[280:361,160:321],hgt2016[0,280:361,160:321],
transform=proj,levels=contourline,linewidths=1,colors="k")
cs3.clabel(contourline[0::2],fontsize=fontNum,colors="k",fmt='%.f')
ax3.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'],alpha=0.5)
ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax3.add_feature(cfeature.STATES.with_scale('50m'))
gl3 = ax3.gridlines(crs=proj,ylocs=np.arange(region[2],region[3]+5,5),
xlocs=np.arange(region[0],region[1]+5,5),draw_labels=False)
ax3.set_title("201601",fontsize=15)
#----第四张图
ax4 = fig.add_axes([0.52,0.32,0.45,0.4],
projection=ccrs.LambertConformal(central_longitude=110.0,
central_latitude=20,standard_parallels=(15,25)))
ax4.set_extent(region,crs=ccrs.PlateCarree())
cs4 = ax4.contour(lons[280:361,160:321],lats[280:361,160:321],hgt2016[1,280:361,160:321],
transform=proj,levels=contourline,linewidths=1,colors="k")
cs4.clabel(contourline[0::2],fontsize=fontNum,colors="k",fmt='%.f')
ax4.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'],alpha=0.5)
ax4.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax4.add_feature(cfeature.STATES.with_scale('50m'))
gl4 = ax4.gridlines(crs=proj,ylocs=np.arange(region[2],region[3]+5,5),
xlocs=np.arange(region[0],region[1]+5,5),draw_labels=False)
ax4.set_title("201602",fontsize=15)
#----第五张图
ax5 = fig.add_axes([0.06,0.02,0.45,0.4],
projection=ccrs.LambertConformal(central_longitude=110.0,
central_latitude=20,standard_parallels=(15,25)))
ax5.set_extent(region,crs=ccrs.PlateCarree())
cs5 = ax5.contour(lons[280:361,160:321],lats[280:361,160:321],hgt2016[2,280:361,160:321],
transform=proj,levels=contourline,linewidths=1,colors="k")
cs5.clabel(contourline[0::2],fontsize=fontNum,colors="k",fmt='%.f')
ax5.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'],alpha=0.5)
ax5.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax5.add_feature(cfeature.STATES.with_scale('50m'))
gl5 = ax5.gridlines(crs=proj,ylocs=np.arange(region[2],region[3]+5,5),
xlocs=np.arange(region[0],region[1]+5,5),draw_labels=False)
ax5.set_title("201603",fontsize=15)
#----第六张图
ax6 = fig.add_axes([0.52,0.02,0.45,0.4],
projection=ccrs.LambertConformal(central_longitude=110.0,
central_latitude=20,standard_parallels=(15,25)))
ax6.set_extent(region,crs=ccrs.PlateCarree())
cs6 = ax6.contour(lons[280:361,160:321],lats[280:361,160:321],hgt2016[3,280:361,160:321],
transform=proj,levels=contourline,linewidths=1,colors="k")
cs6.clabel(contourline[0::2],fontsize=fontNum,colors="k",fmt='%.f')
ax6.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'],alpha=0.5)
ax6.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax6.add_feature(cfeature.STATES.with_scale('50m'))
gl6 = ax6.gridlines(crs=proj,ylocs=np.arange(region[2],region[3]+5,5),
xlocs=np.arange(region[0],region[1]+5,5),draw_labels=False)
ax6.set_title("201604",fontsize=15)
# 加载中国地图
#china = shpreader.Reader('G:/中国地图/maps/bou2_4l.dbf').geometries()
#绘制中国国界省界九段线等等
#ax.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
#---保存
plt.savefig("./高度场Lambert.svg")
plt.savefig("./高度场Lambert.png")
print('-------程序结束-----------')
输出的svg图,几个图会出现地图相互覆盖,
输出png就是正常的。
|
|