Extract temporal data at a special location from grid dataΒΆ
DimArray has interpn
function used to extract a data value with special dimension values. The example
below read the array data from a ARL data file and extract temporal data at a special location.
Extract temporal data
#Set station longitude/latitude
lon = 116.4
lat = 40.3
#Add data file
f = addfile('D:/Temp/arl/gdas1.mar15.w5')
#Read data array
pblh = f['PBLH'][:,:,:]
u = f['U10M'][:,:,:]
v = f['V10M'][:,:,:]
ws = sqrt(u * u + v * v)
#Get time dimension number
tn = pblh.dimlen(0)
#Output data file
ofn = 'D:/Temp/arl/st_test.csv'
of = open(ofn, 'w')
of.write('time,pblh,ws\n')
#Loop
for i in range(tn):
t = f.gettime(i)
#Interpolate to station
pblh_st = pblh[i,:,:].interpn([lat, lon])
ws_st = ws[i,:,:].interpn([lat, lon])
line = '%s,%.2f,%.2f' % (t.strftime('%Y-%m-%d_%H'), pblh_st, ws_st)
print line
of.write(line + '\n')
of.close()
Read extracted data and plot
fn = 'D:/Temp/arl/st_test.csv'
table = readtable(fn, delimiter=',', format='%{yyyy-MM-dd_HH}D%2f')
t = table['time']
pblh = table['pblh']
ws = table['ws']
#Plot
ax1 = axes()
yaxis(ax1, color='r')
plot(t, pblh, 'r-o')
xaxis(axistype='time')
ylabel('Height (m)')
xlim(t[0], t[-1])
title('PBLH and wind speed at location (lon: 116.4, lat: 40.3)')
ax2 = twinx(ax1)
yaxis(ax2, color='b')
plot(t, ws, 'b-^')
ylabel('Wind speed (m/s)')
xlim(t[0], t[-1])