"""
注意！！！
1.列表乘法
a[[]]*3 得到的a中三个列表的地址相同（如果我没有理解错的话）（引用类型会导致意想不到的错误）
2.一些特殊字符的神奇（format函数有更强大的功能）
print('\')得到
print('%%%d'%3) 得到 %3
%-10s 十个字符宽度，左对齐（默认右对齐） %+3d 可以在正数前面加上+号
进制转化：%u 十进制 %0 八进制 %x/%X 十六进制
#的用法：'%#o'%98=='0o142'
0的用法：'%05d'%98=='00098'
%10.2s——10位占位符，截取两位字符串
'%c'%98=='b' #98是ascll码值
3.python中的关键字
False class finally is return None continue for lambda try True
def from nonlocal while and del global not with as elif if or yield
assert else import pass break except in raise
4.理解性问题
x1,x2=[0]*2 那么x1是： 0（int）
"""
from math import *#math库
'''
math 的函数：角度弧度转化：radians(60).degrees(3.14)
指数：exp(5)(e的5次方)
'''
import matplotlib.pyplot as plt #画图库
'''
使matplotlib作图显示中文的动态设置方法（没试过）
from pylab import *
mpl.rcParams['font.sans-serif]=['KaiTi]
'''
import numpy as np #强大的库
from mpl_toolkits.basemap import Basemap #地图绘制库
import warnings
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use('Agg') #这是个啥
import urllib.request as req#从网络上获取文件所需的库
from matplotlib.animation import FuncAnimation #载入动画库
from numpy.linalg import det#用来求行列式的值,linalg是支持线性代数的子包
help(math.sin)#返回对象的帮助信息
dir(math.sin)#返回对象的属性
math.sin.__doc__#返回对象的文档字符串
N_list = [int(i) for i in input().split(',')]
for (vin, rin, rout) in [(10, 10, 5), (10, 10, 10), (20, 20, 5), (20, 20, 10)]:#集体赋值
y.reverse()
t.sort(reverse=true)#排序从大到小，括号里为空则从小到大
round(a,1)#表示a四舍六入五成双，保留一位小数，若没有1则保留整数，-1为十位……
c=complex(a,b)#a/b为c的实/虚部
xstring = xstring.replace('年', '')
xstring = unicode(xstring.encode('cp936').decode('cp936'))#汉字用Unicode编码
a.strip('sda#')#去字符串末尾的一些字符
zip(a,b)#可将两个一维列表打包为二维元组构成的列表（测试失败）
eval(f)#f是一个表达式的字符串形式
a=[1,1,2,11]
a.count(1)==2
a.sort()#无返回值
a[:][0]#这样不行 得不到第一列的数值
outfile=open('file','w')
outfile.write(data)
outfile.close()
'''
sin(radians(x[i]))
y.append()
for i in range(len(a))
y = [f(radians(e)) for e in x]
'''
array类型的一些操作
arange(0,360,1)#开始值，终值（不包括），步长
linspace(0,360,3601)#开始值，终值（包括），元素个数
a=np.array([[1, 2, 3],[4, 5, 6]]).reshape(3,2)#array类型可以加减乘除乘方运算（向量运算）
a=np.mat([[1, 2, 3],[4, 5, 6]])#表示矩阵（mat）
np.hstack()#还没弄明白
np.vstack()#同上
a.sum()
a.max()
c=a.copy()#和c=a不同
a3=np.delete(a2,1,axis=1)#删除二维数组a2中的第二列（axis=1表示列）a2的值不会改变，若无axis则变为一维数组
np.shape(a)#
np.zeros((2,3))=[[0., 0., 0.],[0., 0., 0.]]#可以设置dtype（数据类型）
np.vectorize(H)#向量化(自定义向量化函数可能用到)
ones(3)=[ 1. 1. 1.]
eye(N,M=None, k=0, dtype=<type 'float'>)#第k条对角线全1，其它元素为0的矩阵
identity(3) #结果是3阶单位方阵
a.ndim#数组维数
a.shape#每一维大小
a.size#数组的元素数
a.dtype#元素类型
a.itemsize#每个元素所占的字节数
a[1,1:2]#切片操作，与list注意比较
p=np.polyfit(list1,list2,deg)#用点拟合函数：deg是维数（如3），得到的是一个deg维函数的系数列表如[3,4](表示3*x+4)
f=np.poly1d(p)#将上面得到的列表生成一个函数
画图操作
plt.plot(x,y,'b.')#x,y可以为列表或者array
plt.bar(x,y,width=0.8,color='#ff7f50')#注意rgb与十六进制转换
plt.xlabel('x') #设置x轴文字标记
plt.ylabel('sin/cos') #设置y轴文字标记
plt.legend(['sin(x)', 'cos(x)'], loc='upper center') #图例
"""
x,=plt.plot(xlist,ylist) # , 必不可少（多个返回值）//据说此题还有label的实现方法
plt.legend([x],['sad'],loc="upper left")
"""
plt.grid('on') #显示网格
plt.title("%dnian"%x)
plt.axis([a,b,c,d])
plt.text(x,y,'sda') #纯文本注释
xticks(x, xlabels, rotation = 40)#类似可得yticks
fig=plt.figure(figsize=[12,6])#生成画布
plt.show()
plt.savefig("step2/stu_img/student.png")#保存图表的操作
ani.save('earthquake.mp4', writer='ffmpeg', fps=20)#动画保存
fig= plt.figure(figsize = [12, 6], dpi = 400)#调整清晰度
"""画地图
a="http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/"
b=a+"4.5_month.csv"
d=req.urlopen(b) #个人理解： 生成一个界面
datas=d.read()
datas=str(datas,encoding='ascii',errors='ignore') #有些迷糊
lines=datas.split('
')
data=lines[-2:0:-1]
map=Basemap(projection='robin',lon_0=0,resolution='c') #绘制地图 // map = Basemap(projection='ortho', lat_0=40, lon_0=0)
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.drawmapboundary(fill_color=(0.8,0.95,1.0))
map.fillcontinents(color=(1,0.9,0.7),lake_color=(0.8,0.95,1.0),zorder=0)
for i in range(len(data)):
temp=data[i].split(',')
lat.append(temp[1])
lon.append(temp[2])
for i in range(len(data)):
x1=lat[i]
x2=lon[i]
x,y=map(x2,x1) #经纬度转换为地图上坐标（坐标值非常大）
map.scatter(x, y, s = 40, marker='o', color='r',edgecolor='r')#画点，其实可以整体画一次
plt.text(x, y, u'北京', color = 'b', fontsize=10)#标记汉字’北京‘(u表示中文字符，r表示去掉反斜杠的那些机制，b表示bytes对象，)
#f表示支持{}内python表达式如： print(f'{name} done in {time.time() - t0:.2f} s')
"""
"""自定义向量化函数
一 自动向量化
1.无if检测时相对简单
2.有if检测时
def H(x):
return (0 if x<0 else 1)
H_vec = vectorize(H)
x = arange(-3, 3)
y = H_vec(x)
plot(x, y, '-')
axis([-4,4,-0.1, 1.1])
show()
二 手动向量化
（1）where代替if
from matplotlib.pyplot import plot, show, axis
from numpy import *
def H(x):
return where(x<0, 0.0, 1.0)
x = arange(-3, 3)
y = H(x)
plot(x, y, '-')
axis([-4,4,-0.1, 1.1])
show()
（2）利用循环
"""
"""动画
(1)
def update(frame):
x, y = map(longitude[frame%N], latitude[frame%N])
map.scatter(x, y, s = 40, marker='o', color='r') # 实心红点
animation = FuncAnimation(fig, update, interval=100) # 更新频率单位是千分之一秒
plt.show()
(2)
scat = map.scatter(0,0,s=0,marker='o',color='r') # 初始点位置设为原点，大小取0（即不可见）
def update(frame):
x, y = map(longitude[frame%N], latitude[frame%N]) # 将经纬度转换为图像中的位置
# 修改地震发生点位置和点大小：
scat.set(sizes=[80], offsets=[x, y])#这个函数是什么意思？
animation = FuncAnimation(fig, update, interval=100) # 更新时间为100个千分之一秒 问题：为什么要ani=。。。不加无效果
(3)
pointlimit = 20 # 图中同时最大允许显示的地震次数，相当于缓存的可见点数
# 经度、纬度、地震强度、点大小的缓冲区：
points = np.zeros(4*pointlimit, dtype=float).reshape(pointlimit, 4)
maxsize = 300 # 点大小的最大值，累加超过后从头开始
sizestep = 30 # 点大小一次增加的值
# 颜色的缓冲区：
colors = [[0,0,1,1], [1,0,0,1]] # 两种颜色，前一种表示5.5级及以下，后一种表示5.5级以上
pointcolors = np.zeros(8*pointlimit, dtype=float).reshape(pointlimit, 2, 4) #edgecolor, facecolor
# 动画更新函数：
def update(frame):
recordindex = frame%N # 计算整体数据集中当前处理行数的下标值，处理完之后从头开始
pointindex = frame%pointlimit # 计算缓存可见点的下标值，满了后从头开始
# 所有点的大小均增加一个步长sizeStep
points[:, 3] = [((size+sizestep)%maxsize if size>0 else 0) for size in points[:, 3]]
# 将整体数据集中当前处理数据放入缓存（points）
points[pointindex, 0], points[pointindex, 1] = map(longitude[recordindex], latitude[recordindex])
points[pointindex, 2] = magnitude[recordindex] # 地震强度
points[pointindex, 3] = 20 # 点大小初始化
pointcolors[pointindex,0,:]=(colors[0] if points[pointindex, 2] < 5.5 else colors[1])# edgecolor
pointcolors[pointindex,1,:]=pointcolors[pointindex, 0, :]*(1,1,1,0.25) # facecolor
# 修改散点图的点大小、位置、颜色：
scat.set(sizes=points[:, 3], offsets=zip(points[:, 0],points[:, 1]))
scat.set(edgecolor=pointcolors[:, 0, :], facecolor=pointcolors[:, 1, :])
(4)
scat = map.scatter(0, 0, s = 0, marker='o', color='r') # 初始点位置设为原点，大小取0（即不可见）
ax.set_xlabel(u"近30天大于4.5级强度的地震（红色大于5.5级）", fontsize=16)
ax.set_title('')
pointlimit = 20 # 图中同时最大允许显示的地震次数，相当于缓存的可见点数
# 经度、纬度、地震强度、点大小的缓冲区：
points = np.zeros(4*pointlimit, dtype=float).reshape(pointlimit, 4)
maxsize = 300 # 点大小的最大值，累加超过后从头开始
sizestep = 30 # 点大小一次增加的值
# 颜色的缓冲区：
colors = [[0,0,1,1], [1,0,0,1]] # 两种颜色，前一种表示5.5级及以下，后一种表示5.5级以上
pointcolors = np.zeros(8*pointlimit, dtype=float).reshape(pointlimit, 2, 4) #edgecolor, facecolor
# 动画更新函数：
def update(frame):
recordindex = frame%N # 计算整体数据集中当前处理行数的下标值，处理完之后从头开始
pointindex = frame%pointlimit # 计算缓存可见点的下标值，满了后从头开始
# 所有点的大小均增加一个步长sizeStep
points[:, 3] = [((size+sizestep)%maxsize if size>0 else 0) for size in points[:, 3]]
# 将整体数据集中当前处理数据放入缓存（points）
# 将经纬度投影为地图上的点（map）
points[pointindex, 0], points[pointindex, 1] = map(longitude[recordindex], latitude[recordindex])
points[pointindex, 2] = magnitude[recordindex] # 地震强度
points[pointindex, 3] = 20 # 点大小初始化
pointcolors[pointindex, 0, :] = (colors[0] if points[pointindex, 2] < 5.5 else colors[1]) # edgecolor
pointcolors[pointindex, 1, :] = pointcolors[pointindex, 0, :]*(1,1,1,0.25) # facecolor
# 修改散点图的点大小、位置、颜色：
scat.set(sizes=points[:, 3], offsets=zip(points[:, 0],points[:, 1]))
scat.set(edgecolor=pointcolors[:, 0, :], facecolor=pointcolors[:, 1, :])
# 修改标题以显示动态时间（当前地震发生时间）：
ax.set_title(times[recordindex]) #也许可以直接用的plt.title(times[recordindex]) ，似乎不需要前面的ax.set_title('')
"""