已解决，改为经验，CartoPy画色斑图的时候无法同时显示地图和数据
import numpy as np
from scipy.interpolate import Rbf
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.pyplot import colormaps as cmap
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
from cartopy.mpl.patch import geos_to_path
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
<font face="微软雅黑">data = pd.read_csv('../rain05.csv')
# 插值
lon = data['Lon']
lat = data['Lat']
rain_data = data['PRE_Time_0808']
# [103,110,24.5,29.3]
olon = np.linspace(103,110,88)
olat = np.linspace(24.5,29.3,88)
olon,olat = np.meshgrid(olon,olat)
# 插值处理
func = Rbf(lon, lat, rain_data,function='linear')
rain_data_new = func(olon, olat)
levels = [0,10,25,50,100,150,200,300,450]
fname = '../Micaps-shapefiles/Province.shp'
fname2 = '../Micaps-shapefiles/Province.shp'
adm1_shapes = ShapelyFeature(Reader(fname).geometries(),ccrs.PlateCarree(), facecolor='White',edgecolor='black')
adm2_shapes = ShapelyFeature(Reader(fname2).geometries(),ccrs.PlateCarree(),facecolor='White',edgecolor='black',alpha=0.1,linewidth=2)
#设定裁剪区域
reader = Reader(fname)
provinces = reader.records()
#Micaps地图中Province.shp中可以使用PAC区分不同省份，同理City.shp和County.shp也可使用该关键字
us_multipoly, = [province.geometry for province in provinces if province.attributes['PAC'] == 520000]
main_us_geom = sorted(us_multipoly.geoms, key=lambda geom: geom.area)[0]
us_path, = geos_to_path(main_us_geom)
#画图部分
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
plate_carre_data_transform = ccrs.PlateCarree()._as_mpl_transform(ax)
axis = ax.contourf(olon,olat,rain_data_new, levels=levels, cmap='viridis',transform=ccrs.PlateCarree(),zorder=5)
for collection in axis.collections:
collection.set_clip_path(us_path, plate_carre_data_transform)#设置显示区域
ax.add_feature(adm2_shapes,zorder=1)
# ax.set_extent([103,110,24.5,29.3])</font>
```python
from netCDF4 import Dataset import matplotlib.pyplot as plt import numpy as np import cartopy.io.shapereader as shpreader from cartopy.mpl.patch import geos_to_path from matplotlib.patches import PathPatch import cartopy.crs as ccrs import cartopy.feature as cfeature import shapely.geometry as sgeom import matplotlib.ticker as mticker from cartopy.mpl.gridliner import LONGITUDE_FORMATTER , LATITUDE_FORMATTER import matplotlib.colors as mc
def myset(ax): plate_carre_data_transform = ccrs.PlateCarree()._as_mpl_transform(ax) line=ax.gridlines(draw_labels=True,alpha=0.5,linestyle=’–’) line.xlabels_top=False line.ylabels_left=False #line.xlines=False#只剩下x轴平行线 line.xlocator=mticker.FixedLocator(np.arange(70,140,10))#手动设置x轴刻度 line.xformatter= LONGITUDE_FORMATTER line.yformatter= LATITUDE_FORMATTER line.xlabel_style={‘size’:10,’color’:’black’} line.ylabel_style={‘size’:10,’color’:’black’} line.yalbel_style={‘color’:’red’,’weight’:’bold’} return
def mycon(ax,level,level0,pre,x,y,cmape): plate_carre_data_transform = ccrs.PlateCarree()._as_mpl_transform(ax) file1=r’D:/baiduyundownload/map/China boundary/province boundary/bou2_4p.shp’ reader1 = shpreader.Reader(file1).geometries() for province in reader1: ax.add_geometries([province],ccrs.PlateCarree(),linewidth=0.4,\ edgecolor=’black’,facecolor=’none’,alpha=0.75) norm= mc.BoundaryNorm(level, 256) quad_set = plt.contourf(xlon,xlat,lad,alpha=0.75,levels=level,norm=norm, \ cmap=cmape,transform=ccrs.PlateCarree())#绘制等值线填充 plt.colorbar(quad_set,orientation=’horizontal’,ticks=level0,\ shrink=0.85,extend=’both’,extendrect=False,pad=0.05)
def mlable(ax,level,level0,pre,x,y,cmape): plate_carre_data_transform = ccrs.PlateCarree()._as_mpl_transform(ax) cs=plt.contour(x,y,pre,alpha=0.75,levels=level,colors=’w’,\ linewidths=0.02,\ transform=ccrs.PlateCarree()) cl=plt.clabel(cs,level,fmt=’%i’,colors=’black’,\ fontsize=12,manual=True) for cp in cl: cp.set_rotation(‘horizontal’) cp.set_bbox(dict(facecolor=’white’, edgecolor=’w’, alpha=0.5))#,\ cp.update_bbox_position_size=5
q=Dataset(r’D:/Desktop/ngc_DOMAIN000.nc’) xlon=q.variables[‘xlon’][:] xlat=q.variables[‘xlat’][:] lad=q.variables[‘landuse’][:]
level=range(0,20);level0=range(1,20) cmap=’Dark2’
ax=plt.axes([0.02,0.02,0.9,0.94],projection=ccrs.PlateCarree()) myset(ax) mycon(ax,level,level0,lad,xlon,xlat,cmap) mlable(ax,level,level0,lad,xlon,xlat,cmap)
plt.show()
本文由 xigrug 创作，采用 知识共享署名4.0 国际许可协议进行许可
本站文章除注明转载/出处外，均为本站原创或翻译，转载前请务必署名
最后编辑时间为:2018-07-07 17:03:56