【气象常用】框选研究区域

效果图:

主要步骤:

参考【气象常用】等值线填充+九段线+南海小地图+掩膜中国区域数据-CSDN博客

本次只是加了一个矩形框

详细代码:着急的直接拖到最后有完整代码

步骤一:导入库包及图片存储路径并设置中文字体为宋体,西文为新罗马(没有的库包要先下好奥)

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import xarray as xr
from cnmaps import get_adm_maps, draw_maps
# import pygrib as pg
# 设置西文字体为新罗马字体,中文为宋体,字号为12
from matplotlib import rcParams
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
          }
rcParams.update(config)
rcParams['axes.unicode_minus']=False

步骤二:读入数据、计算距平及掩膜数据

##########################################################
datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'

# 读入数据
f = xr.open_dataset(datapath + 't2m_cn_jja_1950-2023.nc')
# print(f)
t = f['t2m'].loc[:, 55.:0., 70.:140.]
time = f['time']
lons = t['longitude']
lats = t['latitude']
lon, lat = np.meshgrid(lons, lats)
# 计算时间维平均,及距平
t_mean = np.mean(t, axis=0)
dt = t[-1, :, :] - t_mean
# print(np.max(dt))
# print(np.min(dt))

#掩膜数据
cn = get_adm_maps(level='国', only_polygon=True, record='first')
dt = cn.maskout(lon, lat, dt)

步骤三:绘制主图

##########################################################
# 画图
t_level = np.arange(-6, 7, 1)
fig = plt.figure(figsize=(15, 6))
def plot_cont(data, location, level, extent, show):
    ax = fig.add_axes(location, projection=ccrs.PlateCarree(central_longitude=120))
    c1 = ax.contourf(lon, lat, data, levels=level, 
                      extend='both', cmap='RdBu_r', transform=ccrs.PlateCarree())
    # 绘制矩形框  
    rect = plt.Rectangle((105, 28), 120-105, 37-28,  
                          transform=ccrs.PlateCarree(), edgecolor='k', facecolor='none')  
    ax.add_patch(rect)  
    
    # 添加地理信息
    ax.add_feature(cfeature.OCEAN.with_scale('50m'))  # 海洋填充为蓝色
    ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='w')  # 陆地填充为白色
    # help(ax.add_feature)
    # ax.add_feature(cfeature.RIVERS.with_scale('50m'))  # 添加河流
    # ax.add_feature(cfeature.LAKES.with_scale('50m'))  # 添加湖泊
    # ax.add_feature(cfeature.BORDERS.with_scale('50m'),lw=0.5)  # 添加国界线
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=0.5)  # 添加海岸线
    ax.add_geometries(Reader(shppath).geometries(), ccrs.PlateCarree(),
                      facecolor='none', edgecolor='gray', linewidth=0.7)  #添加中国国界

    # 设置显示范围
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    
    # 设置是否显示经纬网格线及坐标刻度
    if show == True:
        # 画经纬度网格
        gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1.2, 
                          color='k', alpha=0.5, linestyle='--')
        
        gl.top_labels = False  # 关闭顶端的经纬度标签
        gl.right_labels = False  # 关闭右侧的经纬度标签
        gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
        gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
    
        # 设置经纬度网格的间隔及坐标标签
        gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+10, 10))
        gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+10, 10))
        ax.set_xticks(list(range(extent[0], extent[1]+10, 10)), crs=ccrs.PlateCarree())
        ax.set_yticks(list(range(extent[2], extent[3]+10, 10)), crs=ccrs.PlateCarree())
        
        lon_formatter = LongitudeFormatter(zero_direction_label=True)  # zero_direction_label用来设置经度的0度加不加E和W
        lat_formatter = LatitudeFormatter()
        ax.xaxis.set_major_formatter(lon_formatter)
        ax.yaxis.set_major_formatter(lat_formatter)
        
    return c1

CN = plot_cont(dt, [0.1, 0.1, 1, 1], t_level, [70, 140, 15, 55], show=True)
NH = plot_cont(dt, [0.68, 0.1, 0.4, 0.4], t_level, [105, 123, 2, 23], show=False)

步骤四:添加色标

###############################################################################
# 色标设置
rect = [0.3, 0.015, 0.6, 0.04]
cbar_ax = fig.add_axes(rect)

cb = fig.colorbar(CN, drawedges=True, ticks=t_level, 
                  cax=cbar_ax, orientation='horizontal')
cb.set_label('℃', fontsize=12)
cb.ax.tick_params(length=0)

步骤五:保存图片

###############################################################################
# 存图
plt.savefig(figpath+'212 框选研究区域', dpi=600, bbox_inches = 'tight')
plt.show()

完整代码在这里:

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import xarray as xr
from cnmaps import get_adm_maps, draw_maps
# import pygrib as pg
# 设置西文字体为新罗马字体,中文为宋体,字号为12
from matplotlib import rcParams
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
          }
rcParams.update(config)
rcParams['axes.unicode_minus']=False
##########################################################
datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'

# 读入数据
f = xr.open_dataset(datapath + 't2m_cn_jja_1950-2023.nc')
# print(f)
t = f['t2m'].loc[:, 55.:0., 70.:140.]
time = f['time']
lons = t['longitude']
lats = t['latitude']
lon, lat = np.meshgrid(lons, lats)
# 计算时间维平均,及距平
t_mean = np.mean(t, axis=0)
dt = t[-1, :, :] - t_mean
# print(np.max(dt))
# print(np.min(dt))

#掩膜数据
cn = get_adm_maps(level='国', only_polygon=True, record='first')
dt = cn.maskout(lon, lat, dt)
##########################################################
# 画图
t_level = np.arange(-6, 7, 1)
fig = plt.figure(figsize=(15, 6))
def plot_cont(data, location, level, extent, show):
    ax = fig.add_axes(location, projection=ccrs.PlateCarree(central_longitude=120))
    c1 = ax.contourf(lon, lat, data, levels=level, 
                      extend='both', cmap='RdBu_r', transform=ccrs.PlateCarree())
    # 绘制矩形框  
    rect = plt.Rectangle((105, 28), 120-105, 37-28,  
                          transform=ccrs.PlateCarree(), edgecolor='k', facecolor='none')  
    ax.add_patch(rect)  
    
    # 添加地理信息
    ax.add_feature(cfeature.OCEAN.with_scale('50m'))  # 海洋填充为蓝色
    ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='w')  # 陆地填充为白色
    # help(ax.add_feature)
    # ax.add_feature(cfeature.RIVERS.with_scale('50m'))  # 添加河流
    # ax.add_feature(cfeature.LAKES.with_scale('50m'))  # 添加湖泊
    # ax.add_feature(cfeature.BORDERS.with_scale('50m'),lw=0.5)  # 添加国界线
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=0.5)  # 添加海岸线
    ax.add_geometries(Reader(shppath).geometries(), ccrs.PlateCarree(),
                      facecolor='none', edgecolor='gray', linewidth=0.7)  #添加中国国界

    # 设置显示范围
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    
    # 设置是否显示经纬网格线及坐标刻度
    if show == True:
        # 画经纬度网格
        gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1.2, 
                          color='k', alpha=0.5, linestyle='--')
        
        gl.top_labels = False  # 关闭顶端的经纬度标签
        gl.right_labels = False  # 关闭右侧的经纬度标签
        gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
        gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
    
        # 设置经纬度网格的间隔及坐标标签
        gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+10, 10))
        gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+10, 10))
        ax.set_xticks(list(range(extent[0], extent[1]+10, 10)), crs=ccrs.PlateCarree())
        ax.set_yticks(list(range(extent[2], extent[3]+10, 10)), crs=ccrs.PlateCarree())
        
        lon_formatter = LongitudeFormatter(zero_direction_label=True)  # zero_direction_label用来设置经度的0度加不加E和W
        lat_formatter = LatitudeFormatter()
        ax.xaxis.set_major_formatter(lon_formatter)
        ax.yaxis.set_major_formatter(lat_formatter)
        
    return c1

CN = plot_cont(dt, [0.1, 0.1, 1, 1], t_level, [70, 140, 15, 55], show=True)
NH = plot_cont(dt, [0.68, 0.1, 0.4, 0.4], t_level, [105, 123, 2, 23], show=False)
###############################################################################
# 色标设置
rect = [0.3, 0.015, 0.6, 0.04]
cbar_ax = fig.add_axes(rect)

cb = fig.colorbar(CN, drawedges=True, ticks=t_level, 
                  cax=cbar_ax, orientation='horizontal')
cb.set_label('℃', fontsize=12)
cb.ax.tick_params(length=0)

###############################################################################
# 存图
plt.savefig(figpath+'212 框选研究区域', dpi=600, bbox_inches = 'tight')
plt.show()

  • 15
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值