MeteoInfoLab脚本示例：站点数据绘制等值线

tonny_xjy

不错，可以试试

强不强

马一下！！！！

ymy_00000

王老师，这是我的程序：
#Set data folders

st_datadir='H:/ymy/Precipitaion6h/2015/12h/06/'
st=datetime.datetime(2015,6,1,12)
et=datetime.datetime(2015,6,1,12)
fns=[]

#Read station data
while st<=et:
  fn=os.path.join(st_datadir,'SURF_WEA_PRE_6HOUR-'+st.strftime('%Y%m%d%H')+'.txt')
  ncol = numasciicol(fn)
  nrow = numasciirow(fn)
  a = asciiread(fn,shape=(nrow,ncol))
  lon= a[:,6]
  lat = a[:,5]
  rain= a[:,8]

#griddata function - interpolate
  x = arange(75, 135, 0.28125)
  y = arange(10, 55, 0.28125)
  prg= griddata((lon, lat), rain, xi=(x, y), method='idw', radius=0.5)
  fns.append(fn)
  st=st+datetime.timedelta(days=1)

#Plot
proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)
axesm(projinfo=proj, position=[0, 0, 0.9, 1], axison=False, gridlabel=False, frameon=False)
geoshow('cn_province', edgecolor='lightgray')
bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
geoshow('cn_cities', facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
china_layer = geoshow('china', visible=False)
levs = [0, 0.1, 10, 25, 50, 100]#降水量级
cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    (240,40,0),(180,10,0),(120,10,0)]
layer = contourfm(x, y, prg, levs, colors=cols)
masklayer(china_layer, [layer])
colorbar(layer, shrink=0.5, aspect=15)
axism([78, 130, 14, 53])#经纬度
text(95, 53, u'全国降水量实况图', fontname=u'黑体', fontsize=18)
text(95, 51, u'(2015-06-01 12:00 至 2015-06-02 12:00)', fontname=u'黑体', fontsize=16)

#Add south China Sea
sc_layer = bou1_layer.clone()
axesm(position=[0.1,0.05,0.15,0.2], axison=False, frameon=True)
geoshow(sc_layer, facecolor=(0,0,255))
xlim(106, 123)
ylim(2, 23)
#savefig('D:/Temp/test/rain_test.png', 800, 800)

ymy_00000

这是我运行程序时候出现的问题：
>>> run script...
Traceback (most recent call last):
  File "<iostream>", line 35, in <module>
  File "E:\ymy\MetoeInfo\MeteoInfo\pylib\mipylib\plotlib\miplot.py", line 2290, in contourfm
    r = gca.contourf(*args, **kwargs)
  File "E:\ymy\MetoeInfo\MeteoInfo\pylib\mipylib\plotlib\mapaxes.py", line 889, in contourf
    ls = plotutil.getlegendscheme(args, a.min(), a.max(), **kwargs)
AttributeError: 'tuple' object has no attribute 'min'

想请问下王老师，这是什么问题导致的？谢谢

MeteoInfo

ymy_00000 发表于 2019-2-27 10:47
王老师，这是我的程序：
#Set data folders

prg= griddata((lon, lat), rain, xi=(x, y), method='idw', radius=0.5)
改为：
prg= griddata((lon, lat), rain, xi=(x, y), method='idw', radius=0.5)[0]

ymy_00000

MeteoInfo 发表于 2019-2-27 12:11
prg= griddata((lon, lat), rain, xi=(x, y), method='idw', radius=0.5)
改为：
prg= griddata((lon,  ...

谢谢王老师！

yuxingzhan

王老师，请问怎么用IDW_Neighbors方式画等值线图？

杜小卜

王老师，你好，出现面图插值不完整该怎么解决？

杜小卜

杜小卜 发表于 2020-4-27 10:29
王老师，你好，出现面图插值不完整该怎么解决？

问题已解决，是因为输出要插值计算的面图，范围x的设置小的问题。

unique

王老师您好，我在画等值线图，但是画出来效果不如Grads简洁，数据已经经过9点平滑，看起来还是有些啰嗦，线条也没有Grads好看（见附件图），想请问王老师如何解决，在此谢过！还有就是等值线label的时候有没有选项说我必须强制标注，或者选择性标注？谢谢！
使用脚本如下：
#Set data folders
#basedir = 'G:\MeteoInfo_2.2.6\mapdata(1)\'
#mapdir = os.path.join(basedir, 'mapdata')
mapdir = 'G:\MeteoInfo_2.2.6\mapdata(1)\mapdata'
#Read shape files
f = addfile('G:/MeteoInfo_2.2.6/nc/temjjadif.nc', 'r')
lai = f['tem']
data=lai[:,:]
f1 = addfile('G:/MeteoInfo_2.2.6/nc/t2difsm.nc', 'r')
lai1 = f1['a']
data1=lai1[:,:]

bou2_layer = shaperead(os.path.join(mapdir, 'China_land.shp'))
bou1_layer = shaperead(os.path.join(mapdir, 'Other_land.shp'))
river_layer = shaperead(os.path.join(mapdir, 'river.shp'))
sea_layer = shaperead(os.path.join(mapdir, 'sea_line.shp'))
ten_layer = shaperead(os.path.join(mapdir, 'tenline.shp'))
#Plot
proj1 = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)
axesm(projinfo=proj1, axison=True, tickfontsize=12, topaxis=True)
geoshow(bou2_layer, edgecolor='black')
geoshow(bou1_layer, edgecolor='black')
geoshow(ten_layer, edgecolor='black')
geoshow(sea_layer, facecolor=(0,160,255),size=0.5)
#geoshow(river_layer, facecolor=(0,0,255),size=0.5)
#axis([90,125,20,54])
levs=[-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0]
col=[(51,0,102),(51,0,204),(51,51,255),(51,102,255),(102,204,255),(255,255,255),(255,204,0),(255,153,0),(255,102,0),(255,51,51),(180,0,0)]
layer = imshowm(data,levs,colors=col,)
tick=[-1,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1]
colorbar(layer, aspect=20, orientation='horizontal', extendrect=False,ticks=tick)
z1=data1
z2=data1
z1[z1>-0.15] = nan
levs1 = [-0.6,-0.2]
col=['black','black','black']
layer1 = contourm(z1,levs1,colors=col,linestyle=':',linewidth=1)
clabel(layer1,fontsize=8)
axism([90, 125, 20, 54])
savefig("D://t2obs.jpg",dpi=300)