MeteoInfoLab处理中国地面气候V3.0数据集(附中国30米,90米,250米,500米,1000米DEM地形数据)

本文源码作者：王亚强 研究员 中国气象科学研究院

---

中国地面国际交换站气候资料日值数据集(V3.0)可以在“中国气象数据网”上下载，每个月每个要素一个文件（ASCII），具体格式说明见网站的相关文档。以2012年5月温度数据文件为例，各列分别为站号、纬度、经度、高度、年、月、日、平均温度、最低温度、最高温度、平均温度质量控制码、最低温度质量控制码、最高温度质量控制码。

MeteoInfoLab有readtable()函数，很适合读取这种表格式的文本数据，文件中没有列名，需要设置headlines=-1，可以用colnames参数设置各数据列的名称。数据被读入一个表格后，可以用列名来获取某一列的数据

---

1 sql()函数提取各站点2012.5.1气温数据并绘制平均气温的站点分布图。需要注意的是数据中经纬度是度分格式，要转化为度。

#Read table data
fn = 'F:/RMeteoInfo/data/SURF_CLI_CHN_MUL_DAY-TEM-12001-201205.txt'
col_names = ['Stid','Lat','Lon','Alt','Year','Month','Day','T_Ave','T_Max',\
    'T_Min','Code_Ave','Code_Max','Code_Min']
table = readtable(fn, headerlines=-1, format='%s%3f%3i%3f%3i', colnames=col_names)
#Sql - only select rows with 'Day = 1'
table = table.sql('Day = 1')
#Get data from table
stid = table['Stid']
lon = table['Lon']
d = (lon / 100).astype('int')
m = (lon - d * 100) / 60
lon = d + m
lat = table['Lat']
d = (lat / 100).astype('int')
m = (lat - d * 100) / 60
lat = d + m
temp = table['T_Ave'] * 0.1
figure(figsize=[700,400], newfig=False)
ax=axesm(tickfontsize=11,axis=True)
# axesm()
world = shaperead('F:/RMeteoInfo/data/map/country1.shp') 
lchina = shaperead('F:/RMeteoInfo/data/map/bou2_4p.shp')
lchina2 = shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp')
#geoshow(world)
geoshow(lchina,edgecolor='k',size=0.8)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.8)
china_layer = geoshow('china', visible=False)
layer = scatterm(lon, lat, temp, 20, cmap='grads_rainbow')
masklayer(lchina, [layer])
title('Temperature (2012-05-01) ' ,bold=True,fontsize=11)
colorbar(layer,orientation='horizontal',ticklen=0,extendrect=False,shrink=1,pad=0.01,aspect=50,label=u'Temperature(°C)',bold=True,fontsize=11)
yaxis(tickvisible=True,location='left',tickwidth=2,linewidth=2,ticklength=3) #ticklength刻度线长度，tickwidth刻度线宽度，linewidth边框宽度
yaxis(tickvisible=False,location='right',tickwidth=2,linewidth=2,ticklength=4) #分别调试tick的宽度，边框线宽和tick的长度
xaxis(tickvisible=False,location='top',tickwidth=2,linewidth=2,ticklength=4)
xaxis(tickvisible=True,location='bottom',tickwidth=2,linewidth=2,ticklength=3)
xlim(70, 140)
ylim(15, 55)
xticks(arange(70, 140.1, 10),bold=True,fontsize=11)
yticks(arange(15, 55.1, 10),bold=True,fontsize=11)
xlabel('Longitude',bold=True,fontsize=11)
ylabel('Latitude',bold=True,fontsize=11)
ax.scale_bar(0.2,0.35,width=165,linewidth=1,bold=True,fontsize=10,bartype='alternating_bar')
ax.north_arrow(0.2,0.9,width=40,height=40,linewidth=1,bold=True,fontsize=10)
#Add south China Sea
axesm(position=[0.7305,0.238,0.15,0.2], axison=False, frameon=True)
geoshow(lchina,edgecolor='k',size=0.4)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.2)
layer = scatterm(lon, lat, temp, 20, cmap='grads_rainbow')
masklayer(lchina, [layer])
xlim(106, 123)
ylim(2, 23)
savefig('F:/RMeteoInfo/test/plot134.png',dpi=800)

---

2 绘制平均气温的空间分布图

注：idw空间插值时未考虑台湾省缺测问题

#Read table data
fn = 'F:/RMeteoInfo/data/SURF_CLI_CHN_MUL_DAY-TEM-12001-201205.txt'
col_names = ['Stid','Lat','Lon','Alt','Year','Month','Day','T_Ave','T_Max',\
    'T_Min','Code_Ave','Code_Max','Code_Min']
table = readtable(fn, headerlines=-1, format='%s%3f%3i%3f%3i', colnames=col_names)
#Sql - only select rows with 'Day = 1'
table = table.sql('Day = 1')
#Get data from table
stid = table['Stid']
lon = table['Lon']
d = (lon / 100).astype('int')
m = (lon - d * 100) / 60
lon = d + m
lat = table['Lat']
d = (lat / 100).astype('int')
m = (lat - d * 100) / 60
lat = d + m
temp = table['T_Ave'] * 0.1
figure(figsize=[700,400], newfig=False)
ax=axesm(tickfontsize=11,axis=True)
# axesm()
world = shaperead('F:/RMeteoInfo/data/map/country1.shp') 
lchina = shaperead('F:/RMeteoInfo/data/map/bou2_4p.shp')
lchina2 = shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp')
#geoshow(world)
geoshow(lchina,edgecolor='k',size=0.8)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.8)
china_layer = geoshow('china', visible=False)
#To grid data
x = arange(70, 140, 0.5)
y = arange(15, 58, 0.5)
gdata,gx,gy = griddata((lon, lat), temp, xi=(x, y), method='idw')
layer = contourfm(x, y, gdata,cmap='grads_rainbow')
masklayer(lchina, [layer])
title('Temperature (2012-05-01) ' ,bold=True,fontsize=11)
colorbar(layer,orientation='horizontal',ticklen=0,extendrect=False,shrink=1,pad=0.01,aspect=50,label=u'Temperature(°C)',bold=True,fontsize=11)
yaxis(tickvisible=True,location='left',tickwidth=2,linewidth=2,ticklength=3) #ticklength刻度线长度，tickwidth刻度线宽度，linewidth边框宽度
yaxis(tickvisible=False,location='right',tickwidth=2,linewidth=2,ticklength=4) #分别调试tick的宽度，边框线宽和tick的长度
xaxis(tickvisible=False,location='top',tickwidth=2,linewidth=2,ticklength=4)
xaxis(tickvisible=True,location='bottom',tickwidth=2,linewidth=2,ticklength=3)
xlim(70, 140)
ylim(15, 55)
xticks(arange(70, 140.1, 10),bold=True,fontsize=11)
yticks(arange(15, 55.1, 10),bold=True,fontsize=11)
xlabel('Longitude',bold=True,fontsize=11)
ylabel('Latitude',bold=True,fontsize=11)
ax.scale_bar(0.2,0.35,width=165,linewidth=1,bold=True,fontsize=10,bartype='alternating_bar')
ax.north_arrow(0.2,0.9,width=40,height=40,linewidth=1,bold=True,fontsize=10)
#Add south China Sea
axesm(position=[0.7305,0.238,0.15,0.2], axison=False, frameon=True)
geoshow(lchina,edgecolor='k',size=0.4)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.2)
layer = contourfm(x, y, gdata,cmap='grads_rainbow')
masklayer(lchina, [layer])
xlim(106, 123)
ylim(2, 23)
savefig('F:/RMeteoInfo/test/plot134.2.png',dpi=800)

---

3 sql()函数提取54511站月数据并绘制气温变化图

#Read table data
fn = 'F:/RMeteoInfo/data/SURF_CLI_CHN_MUL_DAY-TEM-12001-201205.txt'
col_names = ['Stid','Lat','Lon','Alt','Year','Month','Day','T_Ave','T_Max',\
    'T_Min','Code_Ave','Code_Max','Code_Min']
table = readtable(fn, headerlines=-1, format='%s%3f%3i%3f%3i', colnames=col_names)
#Sql - only select rows with 'Stid = "54511"'
table = table.sql('Stid = "54511"')
#Get data from table
day = table['Day']
t_ave = table['T_Ave'] * 0.1
t_min = table['T_Min'] * 0.1
t_max = table['T_Max'] * 0.1
#Plot
plot(day, t_ave, 'r-o', label=u'平均温度')
plot(day, t_min, 'g-^', label=u'最低温度')
plot(day, t_max, 'b-s', label=u'最高温度')
xlim(1, 31)
ylim(10, 40)
legend(loc='upper center', fontname=u'宋体', orientation='horizontal', frameon=False)
xlabel(u'日', fontname=u'宋体')
ylabel(u'温度（摄氏度）', fontname=u'宋体')
xaxis(minortick=True)
yaxis(minortick=True)
title(u'站点54511日温度变化（2012年5月）', fontname=u'黑体')
savefig('F:/RMeteoInfo/test/plot134.1.png',dpi=800)

---

获取本文中国30米，90米，250米，500米，1000米DEM地形数据的途径：气象水文科研猫公众号后台回复：“中国DEM地形数据”，免费获取百度云免费下载链接。

---