MeteoInfoLab处理MICAPS第一类数据(附1980~2020年中国区域30米和1000米土地利用数据)

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

---

这里演示从MICAPS第一类数据中读取站点经纬度及某一个变量的值（比如温度），然后绘制散点图以及插值为格点数据绘制等值线图。打开MICAPS数据的命令是addfile_micaps()，然后就可以从文件变量里读取各个变量的数据。

---

1 用站点数据绘制散点密度图

f = addfile_micaps('F:/RMeteoInfo/data/data3/MICAPS/10101414.000')
data = f['Temperature'][:]
lon = f['Longitude'][:]
lat = f['Latitude'][:]
t = f.gettime(0)
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.3)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.3)
china_layer = geoshow('china', visible=False)
levs = arange(0, 35, 2)
layer = scatterm(lon, lat, data, levs, cmap='grads_rainbow')
masklayer(lchina, [layer])
title(u'Temperature (' + t.strftime('%Y-%m-%d %H:00') + ')',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.2)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.1)
layer = scatterm(lon, lat, data, levs, cmap='grads_rainbow')
masklayer(lchina, [layer])
xlim(106, 123)
ylim(2, 23)
savefig('F:/RMeteoInfo/test/plot135.1.png',dpi=800)

---

2 用站点数据绘制等值线图

先将站点数据IDW插值(反距离权重空间插值)为格点数据。

f = addfile_micaps('F:/RMeteoInfo/data/data3/MICAPS/10101414.000')
data = f['Temperature'][:]
lon = f['Longitude'][:]
lat = f['Latitude'][:]
t = f.gettime(0)
#To grid data
x = arange(70, 140, 0.5)
y = arange(15, 58, 0.5)
gdata,gx,gy = griddata((lon, lat), data, xi=(x, y), method='idw')
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.3)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.3)
china_layer = geoshow('china', visible=False)
levs = arange(0, 35, 2)
layer = contourfm(x, y, gdata, levs, cmap='grads_rainbow')
masklayer(lchina, [layer])
title(u'Temperature (' + t.strftime('%Y-%m-%d %H:00') + ')',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.2)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.1)
layer = contourfm(x, y, gdata, levs, cmap='grads_rainbow')
masklayer(lchina, [layer])
xlim(106, 123)
ylim(2, 23)
savefig('F:/RMeteoInfo/test/plot135.3.png',dpi=800)

---

获取本文1980年、1990年、1995年、2000年、2005年、2010年、2015年、2020年中国区域30米和1000米土地利用数据的途径：气象水文科研猫公众号后台回复：“1980~2020年中国土地利用数据”，免费获取百度云免费下载链接。

---