- 积分
- 3259
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2017-11-10
- 最后登录
- 1970-1-1
|
发表于 2020-10-21 11:14:47
|
显示全部楼层
- import maskout
- import pandas as pd
- import numpy as np
- import matplotlib.pyplot as plt
- from scipy.interpolate import Rbf
- from mpl_toolkits.basemap import Basemap
- plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
- plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
- data = pd.read_excel('D:\climate.xlsx', header=0,
- names=['lon','lat','data'] )
- # 插值
- lon = data['lon']
- lat = data['lat']
- rain_data = data['data']
- olon = np.linspace(97,118,30)
- olat = np.linspace(17,30,30)
- olon,olat = np.meshgrid(olon,olat)
- # 插值处理
- func = Rbf(lon, lat, rain_data,function='linear')
- rain_data_new = func(olon, olat)
- # 画图
- fig = plt.figure(figsize=(16,9))
- plt.rc('font',size=15,weight='bold')
- ax = fig.add_subplot(111)
- m = Basemap(projection='cyl',llcrnrlat=17,llcrnrlon=97,urcrnrlat=30,urcrnrlon=118)
- m.readshapefile('D:/shp_file/four_province','whatevername',color='black')
- x,y = m(olon,olat)
- xx,yy = m(lon,lat)
- levels = np.linspace(0,np.max(rain_data_new),100)
- cf = m.contourf(x,y,rain_data_new, levels=levels)
- cbar = m.colorbar(cf,location='right',format='%d',size=0.3,ticks=np.linspace(0,np.max(rain_data_new),6))
- lon_num = np.arange(97,120,3)
- lon_label = ['97°','99°','102°','105°','108°','111°','114°','117°E']
- lat_num = np.arange(17,32,3)
- lat_label = ['17°','20°','23°','26°','29°N']
- plt.yticks(lat_num,lat_label)
- plt.xticks(lon_num,lon_label)
- plt.title('测试图')
- # 白化
- clip = maskout.shp2clip(cf,ax,'D:/shp_file/four_province',[440000,450000,460000,530000])
- plt.savefig('test.png', bbox_inches='tight',dpi=300)
- plt.show()
复制代码
楼主 请教一下 我这个代码 把特征值从0-9都试了一遍 还是报错local variable 'clip' referenced before assignment 我实在不知道为什么了
这个是maskout的
- import shapefile
- from matplotlib.path import Path
- from matplotlib.patches import PathPatch
- def shp2clip(originfig,ax,shpfile,region):
- sf = shapefile.Reader(shpfile)
- vertices = [] ######这块是已经修改的地方
- codes = [] ######这块是已经修改的地方
- for shape_rec in sf.shapeRecords():
- # if shape_rec.record[3] == region: ####这里需要找到和region匹配的唯一标识符,record[]中必有一项是对应的。
- if shape_rec.record[5] in region: ######这块是已经修改的地方
- 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
复制代码
我的shp文件使用QGIS导出的 我是真的不知道原因在哪里 希望能得到解答 谢谢 |
|