作者： PurpoolPoolObservatory

import matplotlib.pyplot as plt

import matplotlib.dates as mdates

import datetime



years=mdates.YearLocator()

months=mdates.MonthLocator()

yearsFmt=mdates.DateFormatter('%Y')

mag_min=np.zeros(len(mag_array))

fig,ax=plt.subplots(nrows=2,ncols=1,figsize=(6,8))

ax[0].vlines(date_array,mag_min,mag_array,lw=0.1) ### 绘制m-t图

ax[0].set_xlim(datetime.datetime(2009,1,1),datetime.datetime(2018,2,28))

ax[0].set_ylim(0,np.ceil(max(mag_array)))

ax[0].xaxis.set_major_locator(years)  ### x坐标主标记刻度

ax[0].xaxis.set_major_formatter(yearsFmt)  ### x坐标主标记格式

ax[0].xaxis.set_minor_locator(months) ### x坐标次刻度

ax[0].set_xlabel('time')

ax[0].set_ylabel('magnitude')

​

ax[1].bar(date_markers,date_freq,width=35) ### 绘制n-t图，用bar绘制柱状图，这里data_markers是datetime格式的月份，因此width要30天左右

ax[1].set_xlim(datetime.datetime(2009,1,1),datetime.datetime(2018,2,28))

#ax[1].set_ylim(0,np.ceil(max(date_freq)))

ax[1].set_xlabel('time')

ax[1].set_ylabel('monthly freq')

plt.show()

https://github.com/rveciana/geoexamples
import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

def hillshade(array, azimuth, angle_altitude):  # 刻画视觉3D效果

    x, y = gradient(array)

    slope = pi/2. - arctan(sqrt(x*x + y*y))

    aspect = arctan2(-x, y)

    azimuthrad = azimuth*pi / 180.

    altituderad = angle_altitude*pi / 180.

    shaded = sin(altituderad) * sin(slope)\

     + cos(altituderad) * cos(slope)\

     * cos(azimuthrad - aspect)

    return 255*(shaded + 1)/2

fig,ax=plt.subplots()

m=Basemap(projection='merc',llcrnrlat=29,llcrnrlon=102,urcrnrlat=30,urcrnrlon=104,ax=ax)

#ds = gdal.Open('/home/loong1/data/Asia_topo30.grd')

ds = gdal.Open('leshan.grd') #事先准备好grd数据文件

band = ds.GetRasterBand(1)

arr = band.ReadAsArray()#可以从grd文件里读取部分内容，但这里读取了全部

lon,lat=103,29.5

x,y=m(lon,lat)

hs_array = hillshade(arr,315, 45)

m.scatter(x,y,marker='o',c='r')

(x1,y1)=m(102,29)

(x2,y2)=m(104,30)

m.imshow(hs_array,cmap='Greys',origin='upper',extent=[x1,x2,y1,y2]) # extent参数指定图片的坐标范围，叠加到了basemap给定的投影底图上

m.drawparallels(arange(29,30,0.5),labels=[0,1,1,0])

m.drawmeridians(arange(102,104,1),labels=[1,0,0,1])

plt.show()

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple xxxx

自带的cmap本身就有反转版本，一般会在色标名称后面加_r，比如

cm=plt.cm.get_cmap('seismic')

它的反转版本就是

cm=plt.cm.get_cmap('seismic_r')

比如有这样的画图语句

sc=plt.scatter(x,y,c=depth,cmap=cm,...)

cbar=plt.colorbar(sc)

cbar.ax.set_ylabel('depth/km')

参考自/questions/21566610/crop-out-partial-image-using-numpy-or-scipy

判断某点是否处于某个多边形内或边界上在matlab中可以方便地用inpolygon函数达到，python的某些模块中有类似的功能。

import numpy as np

from matplotlib.path import Path

xc = np.array([219.5,284.8,340.8,363.5,342.2,308.8,236.8,214.2])

yc = np.array([284.8,220.8,203.5,252.8,328.8,386.2,382.2,328.8])

xycrop=np.vstack((xc,yc)).T #这里的vstack就是垂直叠加组合形成一个新的数组，T是转置

rnd=np.random.randn(50)*50.0+250.0

rnd=rnd.reshape(25,2)

pth=Path(xycrop,closed=False)

mask=pth.contains_points(rnd) # mask里包含了对应点是否包含在多边形中的逻辑数值

可以用rnd[mask,:]获得包含在多边形内的点，用rnd[~mask,:]获得在外面的点。

perl中是否有matlab中类似于zeros那样的可快速初始化数组的语句？答案是有的，可以这么写：

@a=(0) x 50;

即可生成一个1x50的数据为0的数组。

from dateutil import parser

a=parser.parse('20100101123145')

 

>>a=datetime(2010,1,1,12,31,45)

from dateutil import parser 

a=parser.parse('20100101')

如果想获得这个时间的字符串表达，则可用

a.strftime('%Y%m%d')

https://s tackoverflow.com

/questions/6618515/sorting-list-based-on-values-from-another-list

比如有如下两个list：

people =['Jim','Pam','Micheal','Dwight']

ages =[27,25,4,9]

按ages排序

方法一：

[x for _,x in sorted(zip(ages,people))]

方法二：

import numpy as np

people=np.array(people)

ages=np.array(ages)

indx=ages.argsort()

sortedpeople=people[indx]

方法三：

sorted_index=sorted(range(len(age)),key=lambda x:people[x])

[people[i] for i in sorted_index]



顺便鄙视一下百度和墙。