pandas日期及添加列

已经有了pandas数据框格式的地震目录,需要获取日期的datetime数值,并存储到新的一列dt中,可以这么操作
cata['dt']=cata.apply(lambda x: datetime(x['yr'],x['mn'],x['dy']),axis=1)
这个是按行操作的,很慢。
或者用pandas的to_datetime方法:
pd.to_datetime(df['year'].astype(str) + '-'+ df['month'].astype(str)+ '-1')

而批量操作不行,比如
datetime(x['yr'],x['mn'],x['dy'])
因为datetime不支持pandas的批量操作。

pandas生成日期序列

pandas的date_range函数可以方便的生成日期序列

from datetime import datetime​

import pandas as pd

t_series=pd.date_range(start=datetime(2000,1,1),end=datetime(2000,12,31),freq='D')

这里生成从start到end按天的日期序列,如果freq是M,则是每个自然月的月底​

需要注意的是,现在的t_series并不是datetime类型的序列,而是​timestamp类型的,因此在使用的时候还要注意转换,采用map和匿名函数可以进行方便的转换

t_series_new=t_series.map(lambda t:t.date())​

此时的t_series_new就是date类型的了,还需要注意datetime和date类型是不可以进行比较的,要转成一致才可以。​

matplotlib绘制日期序列图像

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()

如何用matplotlib画dem图像

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()