"""                  
注意！！！
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('\n')
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('')

"""