日志
python 绘制南海和中国地图
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 12 16:11:15 2021
@author: cl
"""
import numpy as np
import pandas as pd
from scipy.interpolate import Rbf # 径向基函数 ： 将站点信息插到格点上 用于绘制等值线
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib as mpl
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from matplotlib.font_manager import FontProperties
import maskout
import os
from matplotlib import font_manager as fm, rcParams
proj = ccrs.PlateCarree(central_longitude=65)
data_2007 = pd.read_table(r'd:/out_data/rt_ano_07_sta.csv',index_col=0,sep=',')
data_1981 = pd.read_table(r'd:/out_data/rt_ano_81_sta.csv',index_col=0,sep=',')
lon = data_1981['lon']
lat = data_1981['lat']
olon = np.linspace(95, 130, 700)
olat = np.linspace(15, 40, 500)
olon, olat = np.meshgrid(olon, olat)
rain_2007 = data_2007['temp']
rain_1981 = data_1981['temp']
func1 = Rbf(lon, lat, rain_2007, function='linear')
func2 = Rbf(lon, lat, rain_1981, function='linear')
rain_2007_new = func1(olon,olat) #插值
rain_1981_new = func2(olon,olat)
# rain_2007_new[rain_2007_new <-50 ] = 0
# rain_1981_new[rain_1981_new <-50 ] = 0
fig = plt.figure(figsize = (10,6),facecolor = '#666666',edgecolor = 'Blue',frameon = False)
ax1 = fig.add_subplot(121,projection=ccrs.PlateCarree())
ax2 = fig.add_subplot(122,projection=ccrs.PlateCarree())
y=[22.5,22.5,32.5,32.5,22.5]
x=[100,110,110,100,100]
levels = np.arange(-4,4,0.2)
cnorm = colors.DivergingNorm(vmin=-4, vcenter=0, vmax=4)
r1 = ax1.contourf(olon,olat,rain_2007_new,levels,transform = ccrs.PlateCarree(),cmap="RdBu_r",norm = cnorm, extend='both')
ax1.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.set_xticks(np.arange(95,130,10),crs = ccrs.PlateCarree()) # x轴
ax1.set_yticks(np.arange(15,40,10),crs = ccrs.PlateCarree()) # y轴
ax1.plot(x,y,lw=1,ls="-",c='black',transform = ccrs.PlateCarree())
fpath = os.path.join(rcParams["datapath"], "fonts/ttf/cmr10.ttf")
prop = fm.FontProperties(fname=fpath)
fname = os.path.split(fpath)[1]
ax1.set_title('(a)2007'.format(fname),loc='left',fontsize=12)
china = shpreader.Reader('D:/shp/map/gadm36_CHN_shp/gadm36_CHN_0.dbf').geometries()
shen=shpreader.Reader('D:/shp/Province_9/Province_9.shp').geometries()
ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
ax1.add_geometries(shen, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
clip = maskout.shp2clip(r1,ax1,'d:/MeteoInfo/map/country.shp','CH') #haerbin/shijie.shp shp/黑龙江省/市界.shp
plt.tick_params(labelsize=12)
labels = ax1.get_xticklabels() + ax1.get_yticklabels()
ax1.spines['bottom'].set_linewidth('1.5')
[label.set_fontname('Times New Roman') for label in labels]
ax = fig.add_axes([0.392, 0.29, 0.1, 0.13],projection=ccrs.PlateCarree())
ax.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
# ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
# 添加南海子图
ax= fig.add_axes([0.392, 0.29, 0.1, 0.13],projection=ccrs.PlateCarree())
ax.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
# ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('D:/shp/map/gadm36_CHN_shp/gadm36_CHN_0.dbf').geometries()
shen=shpreader.Reader('D:/shp/Province_9/Province_9.shp').geometries()
ax.add_geometries(shen, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
ax.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
olon1 = np.linspace(105, 125, 700)
olat1 = np.linspace(0, 25, 500)
olon1, olat1 = np.meshgrid(olon1, olat1)
rain_2007 = data_2007['temp']
rain_1981 = data_1981['temp']
# func1 = Rbf(lon, lat, rain_2007, function='linear')
func2 = Rbf(lon, lat, rain_2007, function='linear')
# rain_2007_new = func1(olon,olat) #插值
rain_2007_new1 = func2(olon1,olat1)
china = shpreader.Reader('D:/shp/map/gadm36_CHN_shp/gadm36_CHN_0.dbf').geometries()
shen=shpreader.Reader('D:/shp/Province_9/Province_9.shp').geometries()
r4=ax.contourf(olon1,olat1,rain_2007_new1,levels,transform = ccrs.PlateCarree(),cmap="RdBu_r",norm = cnorm, extend='both')
clip2 = maskout.shp2clip(r4,ax,'d:/MeteoInfo/map/country.shp','CH')
y1=[27.5,27.5,35,35,27.5]
x1=[105,120,120,105,105]
levels = np.arange(-4,4,0.2)
cnorm = colors.DivergingNorm(vmin=-4, vcenter=0, vmax=4)
r2 = ax2.contourf(olon,olat,rain_1981_new,levels,transform = ccrs.PlateCarree(),cmap="RdBu_r",norm = cnorm, extend='both')
# r2 = ax2.contourf(olon,olat,rain_1981_new,clevs,transform = ccrs.PlateCarree(),cmap=my_cmap,norm = norm)
ax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax2.yaxis.set_major_formatter(LatitudeFormatter())
ax2.set_xticks(np.arange(95,130,10),crs = ccrs.PlateCarree()) # x轴
ax2.set_yticks(np.arange(15,40,10),crs = ccrs.PlateCarree()) # y轴
ax2.plot(x1,y1,lw=1,ls="-",c='black',transform = ccrs.PlateCarree())
ax2.set_title('(b)1981',loc='left',fontsize=12)
china = shpreader.Reader('D:/shp/map/gadm36_CHN_shp/gadm36_CHN_0.dbf').geometries()
shen=shpreader.Reader('D:/shp/Province_9/Province_9.shp').geometries()
ax2.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
ax2.add_geometries(shen, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
clip = maskout.shp2clip(r2,ax2,'d:/MeteoInfo/map/country.shp','CH') #haerbin/shijie.shp shp/黑龙江省/市界.shp
plt.tick_params(labelsize=12)
ax2.spines['bottom'].set_linewidth('1.5')
labels = ax2.get_xticklabels() + ax2.get_yticklabels()
[label.set_fontname('Times New Roman') for label in labels]
ax3 = fig.add_axes([0.814, 0.29, 0.1, 0.13],projection=ccrs.PlateCarree())
ax3.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
# ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader('D:/shp/map/gadm36_CHN_shp/gadm36_CHN_0.dbf').geometries()
shen=shpreader.Reader('D:/shp/Province_9/Province_9.shp').geometries()
ax3.add_geometries(shen, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
ax3.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
olon1 = np.linspace(105, 125, 700)
olat1 = np.linspace(0, 25, 500)
olon1, olat1 = np.meshgrid(olon1, olat1)
rain_2007 = data_2007['temp']
rain_1981 = data_1981['temp']
# func1 = Rbf(lon, lat, rain_2007, function='linear')
func2 = Rbf(lon, lat, rain_1981, function='linear')
# rain_2007_new = func1(olon,olat) #插值
rain_1981_new1 = func2(olon1,olat1)
china = shpreader.Reader('D:/shp/map/gadm36_CHN_shp/gadm36_CHN_0.dbf').geometries()
shen=shpreader.Reader('D:/shp/Province_9/Province_9.shp').geometries()
r3=ax3.contourf(olon1,olat1,rain_1981_new1,levels,transform = ccrs.PlateCarree(),cmap="RdBu_r",norm = cnorm, extend='both')
clip2 = maskout.shp2clip(r3,ax3,'d:/MeteoInfo/map/country.shp','CH')
plt.rcParams['font.family'] = 'Times New Roman'
fig.subplots_adjust(right=0.9)
l = 0.92
b = 0.3
w = 0.015
h = 0.4
#对应 l,b,w,h；设置colorbar位置；(左下宽高)
rect = [l,b,w,h]
cbar_ax = fig.add_axes(rect)
plt.colorbar(r2, cax=cbar_ax)
plt.show()
plt.savefig('d:/picture/temp_ano.png',dpi=1000)