# 引用部分

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.crs as ccrs

import cartopy.io.shapereader as shpreader

from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter

import maskout #只显示某个地区

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文

plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

1234567891011121314151617

# 数据准备

df = pd.read_excel('工作簿1.xlsx') #读取Excel

df.columns = ['stid','lon','lat','height','rain']

# df = df[(df['rain']>0)&(df['lat']>35)&(df['lon']>115)&(df['lat']<50)&(df['lon']<130)] #筛选坐标在某个范围内，由于数据量少，直接在EXCEL中做处理

lon = df['lon']

lat = df['lat']

rain = df['rain']

1234567

# 绘图准备

olon = np.linspace (125,131,120) #经纬坐标，0.05°分辨率 118°到126°0.05分辨率是有160个格点

olat = np.linspace (44,47,60) #纬度坐标，0.05°分辨率

olon,olat = np.meshgrid(olon,olat) #生成坐标网格 meshgrid网格化

func = Rbf(lon,lat,rain,function='linear') #插值函数 调用Rbf插值函数中的 cubic 插值法linear

rain_data_new = func(olon,olat) #插值

rain_data_new[rain_data_new <0 ] = 0

1234567

#画布及绘图声明

fig = plt.figure(figsize = (16,9.6),facecolor = '#666666',edgecolor = 'Blue',frameon = False) # 画布

ax = fig.add_subplot(111,projection=ccrs.PlateCarree()) #绘图区

#色彩定制：24小时降水量级色标

clevs = [0.1,10.,25.,50.,100.,250.,500] #自定义颜色列表

cdict = ['#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042'] #自定义颜色列表 '#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042'

my_cmap = colors.ListedColormap(cdict) # 自定义颜色映射 color-map

norm = mpl.colors.BoundaryNorm(clevs,my_cmap.N) # 基于离散区间生成颜色映射索引

# 绘制等值线、等值线填色

cf = ax.contourf(olon,olat,rain_data_new,clevs,transform = ccrs.PlateCarree(),cmap=my_cmap,norm = norm)

# ct = ax.contour(olon,olat,rain_data_new,clevs) # 绘制等值线

# clabel = ax.clabel(ct,fmt = '%i')

position = fig.add_axes([0.82,0.2,0.05,0.2]) #位置【左，下，宽。高】

plt.colorbar(cf,cax=position) # 颜色参照表

position.set_yticklabels((0,10,25,50,100,250,500))

ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))

ax.yaxis.set_major_formatter(LatitudeFormatter())

ax.set_xticks(np.arange(125,131,2),crs = ccrs.PlateCarree()) # x轴

ax.set_yticks(np.arange(44,47,2),crs = ccrs.PlateCarree()) # y轴

# ax.gridlines() #显示背景线

# # # #裁切

# clip = maskout.shp2clip(cf,ax,'haerbin/hrb.shp')

clip = maskout.shp2clip(cf,ax,'haerbin/hrb.shp','hrb') #haerbin/shijie.shp shp/黑龙江省/市界.shp

#从全国实际地图中获取辽宁省的市级边界并加载

# shpname = r'haerbin/hrb.shp' # shp/country1.shp

# adm1_shapes=list(shpreader.Reader(shpname).geometries())

# ax.add_geometries(adm1_shapes[:],ccrs.PlateCarree(),edgecolor='k',facecolor='') #36:72东三省

plt.savefig('111.png')

1234567891011121314151617181920212223242526272829303132

**

maskout.py

**#coding=utf-8

###################################################################################################################################

#####This module enables you to maskout the unneccessary data outside the interest region on a matplotlib-plotted output instance

####################in an effecient way,You can use this script for free ########################################################

#####################################################################################################################################

#####USAGE: INPUT include 'originfig':the matplotlib instance##

# 'ax': the Axes instance

# 'shapefile': the shape file used for generating a basemap A

# 'region':the name of a region of on the basemap A,outside the region the data is to be maskout

# OUTPUT is 'clip' :the the masked-out or clipped matplotlib instance.

import shapefile

from matplotlib.path import Path

from matplotlib.patches import PathPatch

def shp2clip(originfig,ax,shpfile,region):

sf = shapefile.Reader(shpfile,encoding='utf-8')

for shape_rec in sf.shapeRecords():

if shape_rec.record[2] == region: ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。

vertices = []

codes = []

pts = shape_rec.shape.points

prt = list(shape_rec.shape.parts) + [len(pts)]

for i in range(len(prt) - 1):

for j in range(prt[i], prt[i+1]):

vertices.append((pts[j][0], pts[j][1]))

codes += [Path.MOVETO]

codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)

codes += [Path.CLOSEPOLY]

clip = Path(vertices, codes)

clip = PathPatch(clip, transform=ax.transData)

for contour in originfig.collections:

contour.set_clip_path(clip)

return clip

12345678910111213141516171819202122232425262728293031323334

数据(自己造的) 列分别为 ID 、经度 、纬度 、高度(用不到) 、 降水量

还需要一个shp文件，自己百度下载吧。

最终结果如下：哈尔滨地区降水量色斑图

cs