- 积分
- 17505
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-6-7
- 最后登录
- 1970-1-1
|
发表于 2020-5-8 16:41:55
|
显示全部楼层
本帖最后由 洗雨 于 2020-5-14 16:21 编辑
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import requests
import json
import os
import datetime as DT
import time
from datetime import datetime
# 引用部分
import numpy as np
import pandas as pd
from scipy.interpolate import Rbf # 径向基函数 : 将站点信息插到格点上 用于绘制等值线
import test
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib as mpl
from mpl_toolkits.basemap import Basemap
import cmaps
dt1=DT.datetime(2015,5,2,0,0,0)
dt2=DT.datetime(2015,5,2,0,0,0)
utcDt1 = dt1
utcDt2 = dt2
cur0 = pd.date_range(utcDt1,utcDt2,freq='H')
def basemask(cs, ax, map, shpfile):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
# if shape_rec.record[4] in shpfile:
if shape_rec.record[0] >= 0:
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, prt[i+1]):
vertices.append(map(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform = ax.transData)
for contour in cs.collections:
contour.set_clip_path(clip)
print('1234')
return clip
for time1 in cur0:
fn = time1.strftime('%Y%m%d%H')
fid = r'E:\precipitation_test\hour_1\{0}.txt'.format(fn)
#print(fid = input("请输入需要画图的时次和文本格式(r'E:\precipitation_test\hour_1\2015050100.txt'):"))
name = ['station','lon','lat','PRE_1h','TEM']
data = pd.read_csv(fid,sep='\s+',header=None,names=name)
# print(data)
######插值#########
lon = data['lon']
lat = data['lat']
rain_data = data['PRE_1h']
tem_data = data['TEM']
olon = np.linspace(78,100,88)
olat = np.linspace(26,38,88)
olon,olat = np.meshgrid(olon,olat)
######插值处理######
func = Rbf(lon,lat,tem_data,function = 'linear')
tem_data_new = func(olon,olat)
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
#画布及绘图声明
fig = plt.figure(figsize = (16,9.6)) # 画布
# plt.rc('font',size=15,weight='bold')
ax = fig.add_subplot(111) #绘图区
m = Basemap(projection='cyl',llcrnrlat=26,llcrnrlon=78,urcrnrlat=38,urcrnrlon=100)
m.readshapefile('E:/precipitation_test/Tibet_shps/tibet_shp/xizang_all','xizang_all.shp',default_encoding='gbk',drawbounds=True,linewidth=1, color='k')
m.readshapefile('E:/precipitation_test/Tibet_shps/tibet_shp/river_1','river_1.shp',default_encoding='gbk',linewidth=1,color ='b')
m.readshapefile('E:/precipitation_test/Tibet_shps/tibet_shp/river_2','river_2.shp',default_encoding='gbk',linewidth=0.8,color ='b')
m.readshapefile('E:/precipitation_test/Tibet_shps/tibet_shp/river_3','river_3.shp',default_encoding='gbk',linewidth=0.6,color ='b')
m.readshapefile('E:/precipitation_test/Tibet_shps/tibet_shp/xzlake','xzlake.shp',default_encoding='gbk',linewidth=1,color ='b')
x,y = m(olon,olat)
xx,yy = m(lon,lat)
levels = np.linspace(0,np.max(tem_data_new),50)
#色彩定制:24小时降水量级色标
clevs = [-15.,-10.,-7.,-4.,-2.0,0.,2.,4.,7.,10.,15.] #自定义颜色列表
# 绘制等值线、等值线填色
ct = ax.contourf(x,y,tem_data_new,levels=clevs,cmap='seismic')
plt.title('Tibetan Distribution of Temperature')
position = fig.add_axes([0.91,0.15,0.01,0.7]) #位置【左,下,宽,高】
plt.colorbar(ct,cax=position,label='Temperature (℃)') #颜色参照表
# 白化
basemask(ct, ax, m, 'southwest')
print('hahaha')
runfile('E:/precipitation_test/tem_test.py', wdir='E:/precipitation_test')
Traceback (most recent call last):
File "E:\precipitation_test\tem_test.py", line 114, in <module>
basemask(ct, ax, m, 'southwest')
File "E:\precipitation_test\tem_test.py", line 32, in basemask
for shape_rec in sf.shapeRecords():
File "D:\ProgramData\Annaconda3\lib\site-packages\shapefile.py", line 1039, in shapeRecords
for rec in zip(self.shapes(), self.records())])
File "D:\ProgramData\Annaconda3\lib\site-packages\shapefile.py", line 1012, in records
r = self.__record(oid=i)
File "D:\ProgramData\Annaconda3\lib\site-packages\shapefile.py", line 987, in __record
value = u(value, self.encoding, self.encodingErrors)
File "D:\ProgramData\Annaconda3\lib\site-packages\shapefile.py", line 104, in u
return v.decode(encoding, encodingErrors)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xc7 in position 0: invalid continuation byte
raceback (most recent call last):
File "E:\precipitation_test\tem_test.py", line 18, in <module>
import maskout
ValueError: source code string cannot contain null bytes
runfile('E:/precipitation_test/tem_test.py', wdir='E:/precipitation_test')
Traceback (most recent call last):
File "E:\precipitation_test\tem_test.py", line 18, in <module>
import maskout
ModuleNotFoundError: No module named 'maskout'
运行clip = maskout.shp2clip(ct,ax,'E:\\precipitation_test\\Tibet_shps\\bou2_4p',[540000])报错应该怎么解决,我把utf-8改了还是报错 |
|