短期气候预测实习1

xiaoxiao啊

基于jupyter notebook

# 导入所需库

import numpy as np

import pandas as pd

import xarray as xr

import matplotlib.pyplot as plt

import cartopy.crs as ccrs

import cartopy.feature as cfeature

from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

# 更改字体

plt.rcParams['font.sans-serif'] = ['Times New Roman']

plt.rcParams['axes.unicode_minus'] = False

# 处理数据

f = xr.open_dataset('hgt.mon.mean.nc')  # 读取数据集

hgt48_07 = f.hgt.loc[f.time.dt.month.isin([1]), 500, 80:-20].loc['1948':'2007']

hgt08 = f.hgt.loc['2008-01', 500, 80:-20].squeeze()  # 2008年1月数据

x1 = hgt48_07.mean('time')  # 1月（1948-2007）平均

x2 = hgt08 - hgt48_07.mean('time')  # 2008年1月距平

x3 = hgt08 - hgt08.mean('lon')  # 2008年1月纬偏值

# 绘画气候态

def clim_plot():

    fig = plt.figure(figsize=(6, 3), facecolor='w')

    ax = fig.add_axes([0, 0.1, 0.95, 0.9], projection=ccrs.PlateCarree(180))

    ax.set_title(' Jan(1948-2007) Ave 500hPa height', loc='left', fontsize=13)

    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.6, edgecolor='k')

    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())

    ax.set_yticks(np.arange(-20, 100, 20), crs=ccrs.PlateCarree())

    ax.xaxis.set_major_formatter(LongitudeFormatter())

    ax.yaxis.set_major_formatter(LatitudeFormatter())

    ax.tick_params(which='major', width=0.5, length=5)

    ax.set_aspect(1.5)

    ax.spines['geo'].set_linewidth(0.5)

    m = ax.contourf(f.lon, f.lat.loc[80:-20], x1, cmap='jet', levels=np.arange(5000, 6000, 100), transform=ccrs.PlateCarree())

    cax = fig.add_axes([0.1, -0.08, 0.75, 0.06])

    cb = fig.colorbar(m, cax=cax, orientation='horizontal', ticks=np.arange(5100, 5900, 100))

    cb.outline.set_linewidth(0.1)

    cb.ax.tick_params(length=0)

clim_plot()

# 绘画距平场

def ano_plot():

    levels = list(np.arange(-150, 180, 30))

    levels.remove(0)

    fig = plt.figure(figsize=(6, 3), facecolor='w')

    ax = fig.add_axes([0, 0.1, 0.95, 0.9], projection=ccrs.PlateCarree(180))

    ax.spines['geo'].set_linewidth(0.5)

    ax.set_title(' Jan2008 500hPa height Ano', loc='left', fontsize=13)

    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.6, edgecolor='k')

    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())

    ax.set_yticks(np.arange(-20, 100, 20), crs=ccrs.PlateCarree())

    ax.xaxis.set_major_formatter(LongitudeFormatter())

    ax.yaxis.set_major_formatter(LatitudeFormatter())

    ax.tick_params(which='major', width=0.5, length=5)

    m = ax.contourf(f.lon, f.lat.loc[80:-20], x2, cmap='bwr', extend='both', levels=levels, transform=ccrs.PlateCarree())

    ax.set_aspect(1.5)

    cax = fig.add_axes([0.1, -0.05, 0.75, 0.04])

    cb = fig.colorbar(m, cax=cax, orientation='horizontal', ticks=levels)

    cb.outline.set_linewidth(0.1)

    cb.ax.tick_params(length=0)

ano_plot()

# 绘画纬偏场

def zonal_depart_plot():

    levels = list(np.arange(-250, 300, 50))

    levels.remove(0)

    fig = plt.figure(figsize=(6, 3), facecolor='w')

    ax = fig.add_axes([0, 0.1, 0.95, 0.9], projection=ccrs.PlateCarree(180))

    ax.set_title(' Jan2008 500hPa height zonal departure', loc='left', fontsize=13)

    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.6, edgecolor='k')

    ax.set_xticks(np.arange(0, 360, 60), crs=ccrs.PlateCarree())

    ax.set_yticks(np.arange(-20, 100, 20), crs=ccrs.PlateCarree())

    ax.xaxis.set_major_formatter(LongitudeFormatter())

    ax.yaxis.set_major_formatter(LatitudeFormatter())

    ax.tick_params(which='major', width=0.5, length=5)

    ax.spines['geo'].set_linewidth(0.5)

    ax.set_aspect(1.5)

    m = ax.contourf(f.lon, f.lat.loc[80:-20], x3, cmap='bwr', extend='both', levels=levels, transform=ccrs.PlateCarree())

    cax = fig.add_axes([0.1, -0.05, 0.75, 0.04])

    cb = fig.colorbar(m, cax=cax, orientation='horizontal', ticks=levels)

    cb.outline.set_linewidth(0.1)

    cb.ax.tick_params(length=0)

zonal_depart_plot()