Cartopy 教學

本章全部轉載於‘微信公眾號”氣象科研貓“。

主要介紹了Cartopy的13個示例，如下。

1

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

dlon, dlat = 60, 30   #設置步長
xticks = np.arange(0, 360.1, dlon)  #設置繪圖范圍
yticks = np.arange(-90, 90.1, dlat)

fig = plt.figure(figsize=(6,5)) #設置畫布大小
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=180)) #添加圖幅本體
ax.coastlines() #海岸線

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,linewidth=1, linestyle=':', color='k', alpha=0.8) #格網
gl.xlocator = mticker.FixedLocator(xticks)   #格網設置x軸刻度
gl.ylocator = mticker.FixedLocator(yticks)

ax.set_xticks(xticks, crs=ccrs.PlateCarree())  #圖幅設置坐標軸刻度
ax.set_yticks(yticks, crs=ccrs.PlateCarree())

ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))  #設置坐標軸刻度標簽格式
ax.yaxis.set_major_formatter(LatitudeFormatter())
fig_fname = "全球地圖.png"

#plt.savefig(fig_fname, dpi=500, bbox_inches='tight')
plt.savefig('F:/Rpython/lp11/plot2.png',dpi=1200)
plt.show()

2

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
plt.clf()
fig = plt.figure(dpi=200)
# set the projection to PlateCarree
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 1, 1, projection = proj)
ax.set_global()
# set the gridlines to dashed line and the transparency to 0.7
gl = ax.gridlines(ylocs=np.arange(-90,90+30,30),xlocs=np.arange(-180,180+60,60),draw_labels=True,linestyle='--',alpha=0.7)
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
# set background image to 50-natural-earth-1-downsampled.png
ax.stock_img()
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.LAKES)
plt.savefig('F:/Rpython/lp11/plot11.png',dpi=600)
plt.show()

3

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
# Define the interested area
region = [70,140,15,60] # [west,east,south,north]
plt.clf()
fig = plt.figure(dpi=200)
# set the projection to PlateCarree
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 1, 1,projection = proj)
# set the range of the interested area
ax.set_extent(region,crs = ccrs.PlateCarree())
# set the background image to the 1cm:500km raster 
fname = 'F:/Rpython/lp3/HYP_50M_SR_W/HYP_50M_SR_W.tif'
ax.imshow(plt.imread(fname), origin='upper',transform=ccrs.PlateCarree(),extent=[-180, 180, -90, 90])
# set the gridlines to dashed line and the transparency to 0.7
gl = ax.gridlines(ylocs=np.arange(region[2],region[3]+5,5),xlocs=np.arange(region[0],region[1]+10,10),draw_labels=True,linestyle='--',alpha=0.7)
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
railroads = cfeature.NaturalEarthFeature('cultural','railroads','10m',facecolor='None')
urban_areas = cfeature.NaturalEarthFeature('cultural','urban_areas','10m')
# add borders, coastline, rivers, lakes, railroads，and urban areas
ax.add_feature(cfeature.BORDERS.with_scale('50m')) # '50m' for mediate resolution and '10m' for high resolution
ax.add_feature(cfeature.COASTLINE.with_scale('50m')) 
ax.add_feature(cfeature.RIVERS) # low resolution
ax.add_feature(cfeature.LAKES) # low resolution
ax.add_feature(railroads,edgecolor='gray')
ax.add_feature(urban_areas, edgecolor='m')
plt.savefig('F:/Rpython/lp11/plot12.png',dpi=600)
plt.show()

4

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# Read VTEC data from a file
with open('F:/Rpython/lp5/1st-tecmap.dat') as src:
    data = [line for line in src if not line.endswith('LAT/LON1/LON2/DLON/H\n')]
    tec = np.fromstring(''.join(data), dtype=float, sep=' ')
# Generate a geodetic grid
nlats, nlons = 71, 73
lats = np.linspace(-87.5, 87.5, nlats)
lons = np.linspace(-180, 180, nlons)
lons, lats = np.meshgrid(lons, lats)
tec.shape = nlats, nlons
fig = plt.figure(figsize=(8, 3))
# Set projection
ax = plt.axes(projection=ccrs.PlateCarree())
# Plot contour
cp = plt.contourf(lons, lats, tec, transform=ccrs.PlateCarree(), cmap='jet')
ax.coastlines()
# Add a color bar
plt.colorbar(cp)
# Set figure extent & ticks
ax.set_extent([-180, 180, -87.5, 87.5])
ax.set_xticks([-180, -150, -120, -90, -60, -30, 0, 30, 60, 90, 120, 150, 180])
ax.set_yticks([-80, -60, -40, -20,   0,  20,  40,  60,  80])
plt.savefig('F:/Rpython/lp11/plot15.png',dpi=600)
plt.show()

5

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import warnings
import os
from cartopy import config
from pylab import imread
#from cv2  import imread, imshow
def create_map():
    extent = [70, 140, 0, 60]
    shp_path = 'F:/Rpython/lp11/'
    # --創建畫圖空間
    proj = ccrs.PlateCarree()  # 創建坐標系
    fig = plt.figure(figsize=(6, 8), dpi=350)  # 創建頁面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  # 創建子圖
    # --設置地圖屬性
    reader = Reader(shp_path  +  'shp/Province_9.shp')
    provinces = cfeat.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none')
    ax.add_feature(provinces, linewidth=1)
    ax.set_extent(extent, crs=proj)
    # --增加低分辨率地形圖（cartopy自帶）
    #ax.stock_img()  # 增加地形圖
    # --增加高分辨率地形圖（需自行下載）
    #https://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-1/
    #fname = os.path.join(config["repo_data_dir"], 'raster', 'natural_earth', 'NE1_HR_LC.tif')
    fname='F:/Rpython/lp11/shp/NE1_LR_LC.tif'
    ax.imshow(imread(fname), origin='upper', transform=proj, extent=[-180, 180, -90, 90])
    # --設置網格點屬性
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1.2, color='k', alpha=0.5, linestyle='--')
    gl.xlabels_top = False  # 關閉頂端的經緯度標簽
    gl.ylabels_right = 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))
    return ax
if __name__ == '__main__':
    warnings.filterwarnings('ignore')
    ax = create_map()
    plt.savefig('F:/Rpython/lp11/plot17.png',dpi=600)
    plt.show()

6

import warnings
warnings.filterwarnings('ignore')
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure(figsize=(6, 4))
plt.subplots_adjust(left=0.05, right=0.95, top=0.90, bottom=0.05, wspace=0.15, hspace=0.05)
ax = plt.subplot(111)
m = Basemap(resolution='i', projection='merc', llcrnrlat=10.0, urcrnrlat=55.0, llcrnrlon=60., urcrnrlon=140.0)
m.drawcountries(linewidth=0.5)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(10., 55., 10.), labels=[1, 0, 0, 0], color='black',linewidth=0.2,dashes=(None, None))  # draw parallels,dashes=[1,0],
m.drawmeridians(np.arange(60., 140., 10.), labels=[0, 0, 0, 1], color='black',linewidth=0.2,dashes=(None, None))  # draw meridians ,dashes=[1,0]
x, y = m(116.4204, 40.21244)  # Bejing
x2, y2 = m(125.27538, 43.83453)  # Changchun
plt.annotate('Beijing', xy=(x, y), xycoords='data',
              # xytext=(x2, y2), textcoords='offset points',
              xytext=(x2, y2), textcoords='data',
              color='r',arrowprops=dict(arrowstyle="fancy", color='g'))
plt.savefig('F:/Rpython/lp11/plot21.png',dpi=600)
plt.show()

7

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import cmaps
meteo_file = r'F:/Rpython/lp3/hls/ECMWF.nc'
ds = Dataset(meteo_file, mode='r')
# 獲取每個變量的值
lons = ds.variables['longitude'][:]
lats = ds.variables['latitude'][:]
# surface_air_pressure
sp = ds.variables['sp'][:]
sp_units = ds.variables['sp'].units
scale_factor = ds.variables['sp'].scale_factor
add_offset = ds.variables['sp'].add_offset
sp = scale_factor * sp + add_offset
# 2 metre temperature
t2m = ds.variables['t2m'][:]
t2m_units = ds.variables['t2m'].units
scale_factor = ds.variables['t2m'].scale_factor
add_offset = ds.variables['t2m'].add_offset
t2m = scale_factor * t2m + add_offset
# Total column ozone
tco3 = ds.variables['tco3'][:]
tco3_units = ds.variables['tco3'].units
scale_factor = ds.variables['tco3'].scale_factor
add_offset = ds.variables['tco3'].add_offset
tco3 = scale_factor * tco3 + add_offset
# 經緯度平均值
lon_0 = lons.mean()
lat_0 = lats.mean()
# 畫圖大小設置
fig = plt.figure(figsize=(16, 9))
plt.rc('font', size=15, weight='bold')
ax = fig.add_subplot(111)
m = Basemap(lat_0=lat_0, lon_0=lon_0)
lon, lat = np.meshgrid(lons, lats)
xi, yi = m(lon, lat)
# 這里數據時間是UTC 00:00，2018年1月的日平均數據，只展示1月1號的數據
sp_01 = sp[0:1, :, :]
t2m_01 = t2m[0:1, :, :]
tco3_01 = tco3[0:1, :, :]
levels = m.pcolor(xi, yi, np.squeeze(tco3_01), cmap=cmaps.GMT_panoply)
# 添加格網與繪制經緯線
m.drawparallels(np.arange(-90., 91., 20.), labels=[1, 0, 0, 0], fontsize=15)
m.drawmeridians(np.arange(-180., 181., 40.), labels=[0, 0, 0, 1], fontsize=15)
# 添加海岸線,省州邊界以及國家行政邊界
m.drawcoastlines()
m.drawstates()
m.drawcountries()
# 添加colorbar
cbar = m.colorbar(levels, location='bottom', pad="10%")
cbar.set_label(tco3_units, fontsize=15, weight='bold')
# 添加圖的標題
plt.title('Total column ozone')
#----------------------------------------------------------------------------------
plt.savefig('F:/Rpython/lp7/plot22.png',dpi=1200)
plt.show()
ds.close()

8

import numpy as np
import matplotlib.pyplot as plt
import maskout1
from mpl_toolkits.basemap import Basemap
from scipy.interpolate import griddata
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r'F:/Rpython/lp3/guangxi/SimHei/SimHei.ttf')
#創建三個列表用於存放點的數據
x=[]   
y=[]
h=[]
#with open(r'F:/Rpython/lp3/guangxi/廣西站點修正.csv','r') as f:
with open(r'F:/Rpython/lp3/guangxi/廣西海拔高度.csv','r') as f:    
    skip1=f.readline()    
    datas=f.readlines()
for line in datas:
    p=line[:-1]                #去掉每行最后的換行符
    a=p.split(',')             #因為每行每個數據之間用","隔開，去掉逗號。
    x.append(float(a[2]))
    y.append(float(a[3]))
    h.append(float(a[4]))
points=[]
for i in range(len(x)):
    point=[]
    point.append(x[i])
    point.append(y[i])   
    points.append(point)
points=np.array(points)
xi=np.linspace(min(x),max(x),200)
yi=np.linspace(min(y),max(y),200)
xi,yi=np.meshgrid(xi,yi)#網格化
#scipy.interpolate.griddata(points, values, xi, method=‘linear’, fill_value=nan, rescale=Fale)
#zi=griddata(points,h,(xi,yi),method='cubic')
zi=griddata(points,h,(xi,yi),method='linear')
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
# [west,east,south,north]
map = Basemap(
#    llcrnrlon=72,
#    llcrnrlat=0,
#    urcrnrlon=136,
#    urcrnrlat=55,
    llcrnrlon=104,
    llcrnrlat=21,
    urcrnrlon=112.5,
    urcrnrlat=27,    
    resolution = None, 
    projection = 'cyl')
#map.readshapefile('e:\python\map\world0','whatevername',color='black',linewidth=1.2)
#map.readshapefile('e:\python\map\china1','whatevername',color='black',linewidth=0.3)
map.readshapefile(r'F:/Rpython/lp3/guangxi/map/guangxi','whatevername',color='black',linewidth=1.0)
minval,maxval=int(np.amin(h)),int(np.amax(h))
print(minval,maxval)
#c2 = plt.contourf(xi,yi,zi,range(minval,maxval,1),cmap= 'jet')
#c2 = plt.contourf(xi,yi,zi,18,cmap= 'jet')
c2 = plt.contourf(xi,yi,zi,range(0,2000,10),cmap= 'jet')
clip=maskout1.shp2clip(c2,ax,r'F:/Rpython/lp3/guangxi/map/world0',r'F:/Rpython/lp3/guangxi/map/guangxi')
#plt.rcParams['font.family'] = 'Times New Roman'  #設置colorbar的label字體
plt.rcParams['font.family'] = 'DejaVu Sans'   #設置colorbar的label字體
plt.rcParams['font.size'] = 14                #設置colorbar的label字體大小
bar=map.colorbar(c2,label='height / m')
bar.ax.tick_params(labelsize=14)           #設置colorbar刻度字體大小。
# 以下標經緯度
map.drawparallels(np.arange(-90.,90.,2),labels=[1,0,0,0],fontsize=14,linewidth=0.1)
map.drawmeridians(np.arange(-180.,180.,2),labels=[0,0,0,1],fontsize=14,linewidth=0.1)
#u以下寫漢字標題 
plt.title(u'廣西海拔高度',color='red',fontproperties=font,fontsize=14)
#plt.xlabel(u'經度',color='blue',fontproperties=font,fontsize=14)
#plt.ylabel(u'緯度',color='blue',fontproperties=font,fontsize=14)
plt.savefig(r'F:/Rpython/lp3/guangxi/gxfig.png',format='png',dpi=1200)  
plt.show()

9

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.util import add_cyclic_point
#from palettable.colorbrewer.diverging import RdBu_11_r 
import netCDF4 as nc
sns.set_context('talk', font_scale=1.2) # 設置圖形屬性
# read NetCDF file
fn = 'F:/Rpython/lp11/shp/air.sig995.2012.nc'
data = nc.Dataset(fn, 'r') # 默認為讀文件，此處 'r' 可省略
# 讀取相關變量
lat = data.variables['lat'][:].data
lon = data.variables['lon'][:].data
time = data.variables['time'][:].data
air = data.variables['air'][:].data
# 添加數據循環，以防止在0和360時出現白色條狀
# 添加數據循環和不添加數據循環的效果見后文兩張圖
cycle_air, cycle_lon = add_cyclic_point(air, coord=lon)
cycle_LON, cycle_LAT = np.meshgrid(cycle_lon, lat)
projection = ccrs.PlateCarree() # 設置投影
fig, ax = plt.subplots(figsize=(16, 9), subplot_kw=dict(projection=projection))
con = ax.contourf(cycle_LON, cycle_LAT, cycle_air[0, ...],np.arange(220, 321), cmap='Spectral_r')
ax.coastlines(linewidth=1.5) # 添加海岸線
ax.set_xticks(np.arange(-180, 181, 60), crs=projection)
ax.set_yticks(np.arange(-90, 91, 30), crs=projection)
# 設置 ticklabels 格式
lon_formatter = LongitudeFormatter(number_format='.0f',degree_symbol='',dateline_direction_label=True)
lat_formatter = LatitudeFormatter(number_format='.0f',degree_symbol='')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# shrink 控制 colorbar 長度，pad 控制colorbar和圖的距離
cb = fig.colorbar(con, shrink=0.73, pad=0.02)
# 調整 ticklabels
cb.set_ticks(np.arange(220, 321, 20))
cb.set_ticklabels(np.arange(220, 321, 20))
cb.ax.tick_params(direction='in', length=5) # 控制 colorbar tick
# 保存文件, dpi用於設置圖形分辨率, bbox_inches 盡量減小圖形的白色區域
#fig.savefig(F:/Rpython/lp11/plot24.png', dpi=600, bbox_inches='tight')
plt.savefig('F:/Rpython/lp11/plot24.png',dpi=600)
plt.show()

10

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
extent = [70, 140, 0, 60]
Chinese_land_territory_path = 'F:/Rpython/lp11/shp/China_land_territory'
Chinese_10dash_line_path = 'F:/Rpython/lp11/shp/China_10-dash_line'
world_countries_path = 'F:/Rpython/lp11/shp/world_countries'
# 創建坐標系
prj = ccrs.PlateCarree()
fig = plt.figure(figsize=(10, 8), dpi=350)
ax = fig.subplots(1, 1, subplot_kw={'projection': prj})
ax.set_extent(extent, crs=prj)
# 繪制中國陸地領土邊界
    # reader = Reader(shp_path  +  'shp/Province_9.shp')
    # provinces = cfeat.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none')
    # ax.add_feature(provinces, linewidth=1)
Chinese_land_territory2 = Reader('F:/Rpython/lp11/shp/China_land_territory/China_land_territory.shp')
Chinese_land_territory = cfeat.ShapelyFeature(Chinese_land_territory2.geometries(),prj, edgecolor='k',facecolor='none')
ax.add_feature(Chinese_land_territory, linewidth=1)
# 繪制中國十段線
Chinese_10dash_line2 = Reader('F:/Rpython/lp11/shp/China_10-dash_line/China_10-dash_line.shp').geometries()
Chinese_10dash_line = cfeat.ShapelyFeature(Chinese_10dash_line2,prj, edgecolor='r')
ax.add_feature(Chinese_10dash_line, linewidth=2)
# 繪制世界各國領土邊界
world_countries2 = Reader('F:/Rpython/lp11/shp/world_countries/world_countries.shp').geometries()
world_countries = cfeat.ShapelyFeature(world_countries2,prj, edgecolor='k', alpha=0.3,facecolor='none')
ax.add_feature(world_countries, linewidth=0.5)
# 添加自帶地形數據
ax.stock_img()
# 繪制網格點
gl = ax.gridlines(crs=prj, draw_labels=True, linewidth=1.2, color='k',alpha=0.5, linestyle='--')
gl.xlabels_top = False  # 關閉頂端的經緯度標簽
gl.ylabels_right = 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))
plt.savefig('F:/Rpython/lp11/plot25.png',dpi=600)
plt.show()

11

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np 
data = pd.read_csv('F:/Rpython/lp11/shp/1950-2018_all_tornadoes.csv')
fig, ax = plt.subplots(figsize=(12,8))
m = Basemap(projection='mill',llcrnrlat=25,urcrnrlat=50,llcrnrlon=-125,urcrnrlon=-65,rsphere=6371200.,resolution='l',area_thresh=10000)
m.drawcoastlines()
m.drawstates()
parallels = np.arange(25,51,5.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
meridians = np.arange(-125.,-64.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
# 這里要轉換坐標，否則可能無法顯示圖形
#  由於 data 為 DataFrame，所以要直接提取值，因為 basemap 在轉換
# 坐標時只能是列表，元組或標量
x, y = m(data.slon.values, data.slat.values)
hex = m.hexbin(x, y, bins = 150, cmap = plt.get_cmap('Reds'))
fig.colorbar(hex, ax = ax, shrink = 0.6)
plt.savefig('F:/Rpython/lp11/plot28.png',dpi=600)
plt.show()

12

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
# 數據讀取及時間平均處理
ds = xr.open_dataset('F:/Rpython/lp11/shp/EC-Interim_monthly_2018.nc')
temp = (ds['t2m'] - 273.15).mean(dim='time')  #把溫度轉換為℃並對其時間緯求平均
temp.attrs['units'] = 'deg C'  #溫度單位轉換為℃
# 創建畫圖空間
proj = ccrs.PlateCarree()  #創建投影
fig = plt.figure(figsize=(9,6))  #創建頁面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  #子圖
# 設置地圖屬性:加載國界、海岸線、河流、湖泊
ax.add_feature(cfeat.BORDERS.with_scale('50m'), linewidth=0.8, zorder=1)
ax.add_feature(cfeat.COASTLINE.with_scale('50m'), linewidth=0.6, zorder=1)  
ax.add_feature(cfeat.RIVERS.with_scale('50m'), zorder=1)  
ax.add_feature(cfeat.LAKES.with_scale('50m'), zorder=1)  
# 設置網格點屬性
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,linewidth=1.2, color='k', alpha=0.5, linestyle='--')
gl.xlabels_top = False  #關閉頂端標簽
gl.ylabels_right = False  #關閉右側標簽
gl.xformatter = LONGITUDE_FORMATTER  #x軸設為經度格式
gl.yformatter = LATITUDE_FORMATTER  #y軸設為緯度格式
# 設置colorbar
cbar_kwargs = {'orientation': 'horizontal',
   'label': '2m temperature (℃)',
   'shrink': 0.8,
   'ticks': np.arange(-30,30+5,5)}
# 畫圖
levels = np.arange(-30,30+1,1)
temp.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
plt.savefig('F:/Rpython/lp11/plot26.png',dpi=600)
plt.show()

13

import os
import maskout
import numpy as np
import xarray as xr
from scipy.interpolate import Rbf
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from matplotlib.font_manager import FontProperties
SHP = 'F:/Rpython/lp1/wchina/china_shp'
# 數據處理
ds = xr.open_dataset('F:/Rpython/lp1/wchina/r160.nc')
da = ds.PRE.sel(month=slice(6, 8)).mean(dim='month').mean(dim='year')
# 插值
olon = np.linspace(72, 137, 260)
olat = np.linspace(15, 55, 120)
olon, olat = np.meshgrid(olon, olat)
func = Rbf(ds.station_lon, ds.station_lat, da, function='linear')
pre_grid = func(olon, olat)
# 字體
#mpl.rcParams['font.sans-serif'] = ['Times New Roman']
#mpl.rc('font', size=15, weight='normal')
#font = {'family': 'NSimSun','weight': 'normal','size': 15,}
font = FontProperties(fname=r"F:/Rpython/fonts/simsun.ttf")
# 顏色
levels = [0, 0.1, 10, 25, 50, 100, 250, 500]
cmaps = mpl.colors.ListedColormap([ '#FFFFFF', '#A9F18D', '#3DB93D', '#60B8FF', '#0001E2', '#F900FA', '#810045'])
norm = mpl.colors.BoundaryNorm(levels, 7)
# 畫圖
proj = ccrs.LambertConformal(central_latitude=90, central_longitude=105)
fig, ax = plt.subplots(figsize=(15, 15), subplot_kw=dict(projection=proj))
ax.set_extent([80, 130, 13, 55], ccrs.PlateCarree())
ax.gridlines(linestyle='--')
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_geometries(Reader(os.path.join(SHP, 'cnhimap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=1)
cf = ax.contourf(olon,olat,pre_grid,levels=levels, cmap=cmaps,norm=norm,zorder=1,transform=ccrs.PlateCarree())
ax.scatter(ds.station_lon, ds.station_lat,c='k',s=10,marker='o',transform=ccrs.PlateCarree())
cbar = fig.colorbar(cf, ax=ax,orientation="vertical",aspect=25,pad=0.01,shrink=0.7)
cbar.set_label('mm')
cbar.set_ticklabels(levels)
ax.set_title('Annual average summer precipitation')
#  南海
pos1 = ax.get_position()
pos2 = [pos1.x1 - ((pos1.x1 - pos1.x0) / 7), pos1.y0, pos1.width / 7,pos1.height / 5]
proj_n = ccrs.LambertConformal(central_latitude=90, central_longitude=115)
ax_n = fig.add_axes(pos2, projection=proj_n)
ax_n.set_extent([105, 120, 0, 20], ccrs.PlateCarree())
ax_n.add_geometries(Reader(os.path.join(SHP, 'cnhimap.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)
ax_n.add_feature(cfeature.OCEAN.with_scale('50m'))
ax_n.add_feature(cfeature.LAND.with_scale('50m'))
ax_n.add_feature(cfeature.RIVERS.with_scale('50m'))
ax_n.add_feature(cfeature.LAKES.with_scale('50m'))
# 白化
clip = maskout.shp2clip(cf,ax,shpfile=os.path.join(SHP, 'country1.shp'),region='China',proj=proj)
fig.savefig('F:/Rpython/lp11/plot.png', dpi=600, bbox_inches='tight')
plt.show()