探空资料对比
2019-01-01 本文已影响20人
榴莲气象
最后选用METpy画模式和观测的对比（METpy可以兼顾两者画图）
观测的还是选用Wyoming upper air data而不是IDD的，这个垂直层更多
但是IDD的有四个时次的。
查看其它操作Sounding Initialization
image.png
https://bitbucket.org/donaldjmorton/
Intro to Skew-Ts and Sounding Parameters
The sounding class supports the following initialization modes:
- From a University of Wyoming txt file
- From a University of Wyoming Website
- Form an ARM sounding Netcdf file
- From a dictionary with the field names, field values pairs
- Adding individual fields manually
University of Wyoming Sounding Data
Fetch from txt file
The easiest way to get sounding data is to visit the University of Wyoming’s website:
http://weather.uwyo.edu/upperair/sounding.html
To get some sounding data, follow the link and find the date, and location you are interested in, view the data as a text file and just save the file to your system. If you want to get loads of data please be considerate about the way you go about doing this! (Lots of wget requests makes the server administrators unhappy).
You can also pass your own data to SkewT by writing a text file in identical format to the University of Wyoming files, which are constant-width columns. Here’s a sample of the first few lines of one of the bundled examples:
94975 YMHB Hobart Airport Observations at 00Z 02 Jul 2013
PRES HGHT TEMP DWPT RELH MIXR DRCT SKNT THTA THTE THTV
hPa m C C % g/kg deg knot K K K
1004.0 27 12.0 10.2 89 7.84 330 14 284.8 306.7 286.2
1000.0 56 12.4 10.3 87 7.92 325 16 285.6 307.8 286.9
993.0 115 12.8 9.7 81 7.66 311 22 286.5 308.1 287.9
From now on, it’s assumed that the package is installed and the current working directory is the examples one, included in this package.
To read a sounding from a txt file and create a quick plot using the default parameters we only have to do:
from SkewTplus.sounding import sounding
#Load the sounding data
mySounding = sounding("./exampleSounding.txt",fileFormat='txt')
#Do a quick plot
mySounding.quickPlot()
The resulting plot will look like this:
../_images/soundingQuickView.png
Fetch from University Of Wyoming Website
The sounding class supports getting University of Wyoming sound data directly from the UWYO website, we only need to specify the date and the station id.
For example, to initialize the class with with the sounding from “72558 OAX Omaha” station, at April 10th of 2017 OO UTC we simple do:
from SkewTplus.sounding import sounding
#Load the sounding data
mySounding = sounding("20170410:00",fileFormat='web', stationId= "OAX")
#Do a quick plot
mySounding.quickPlot()
ARM Sounding Data
The sounding class also supports initialization from ARM sounding data (Netcdf files). For example:
from SkewTplus.sounding import sounding
#Load the sounding data
mySounding = sounding("./armSoundingExample.cdf",fileFormat='arm')
#Do a quick plot
mySounding.quickPlot()
From a dictionary
The sounding class can be initialized from a dictionary with “field names” , “field values” pairs. The Temperature should be in Celsius degrees and the pressure in hPa.
The next is an example of a dictionary initialization used to plot a sounding from a WRF output file:
from netCDF4 import Dataset
import numpy
from SkewTplus.sounding import sounding
#Load the WRF File
wrfOutputFile = Dataset("wrfOutputExample.nc")
theta = wrfOutputFile.variables["T"][:] + 300 # Potential temperature
# Pressure in hPa
pressure = (wrfOutputFile.variables['P'][:] + wrfOutputFile.variables['PB'][:])
qvapor = wrfOutputFile.variables['QVAPOR'][:]
qvapor = numpy.ma.masked_where(qvapor <0.00002, qvapor)
T0 = 273.15
referencePressure = 100000.0 # [Pa]
epsilon = 0.622 # Rd / Rv
# Temperatures in Celsius
temperature = theta* numpy.power((pressure / referencePressure), 0.2854) - T0
vaporPressure = pressure * qvapor / (epsilon + qvapor)
dewPointTemperature = 243.5 / ((17.67 / numpy.log(vaporPressure * 0.01 / 6.112)) - 1.) #In celsius
dewPointTemperature = numpy.ma.masked_invalid(dewPointTemperature)
# Now we have the pressure, temperature and dew point temperature in the whole domain
# Select one vertical column , t =0 , x=30, y=30
inputData = dict(pressure=pressure[0,:,30,30]/100,
temperature=temperature[0,:,30,30],
dewPointTemperature=dewPointTemperature[0,:,30,30])
mySounding = sounding(inputData)
mySounding.quickPlot()
Adding Fields Manually
The sounding class supports an empty initialization (without any fields). With the addField() method, new fields can be added to the class. With this kind of initialization full control ever the fields added can be obtained. Internally, the class stores the field data values as soundingArrayclasses. This class is a MaskedArray with metadata (long Name,units and missing data value).
To exemplify the use of this initialization, the previous example of the sounding with WRF data coded to use the addField() method:
from netCDF4 import Dataset
import numpy
from SkewTplus.sounding import sounding
#Load the WRF File
wrfOutputFile = Dataset("wrfOutputExample.nc")
theta = wrfOutputFile.variables["T"][:] + 300 # Potential temperature
# Pressure in hPa
pressure = (wrfOutputFile.variables['P'][:] + wrfOutputFile.variables['PB'][:])
qvapor = wrfOutputFile.variables['QVAPOR'][:]
qvapor = numpy.ma.masked_where(qvapor <0.00002, qvapor)
T0 = 273.15
referencePressure = 100000.0 # [Pa]
epsilon = 0.622 # Rd / Rv
# Temperatures in Celsius
temperature = theta* numpy.power((pressure / referencePressure), 0.2854) - T0
vaporPressure = pressure * qvapor / (epsilon + qvapor)
dewPointTemperature = 243.5 / ((17.67 / numpy.log(vaporPressure * 0.01 / 6.112)) - 1.) #In celsius
dewPointTemperature = numpy.ma.masked_invalid(dewPointTemperature)
# Now we have the pressure, temperature and dew point temperature in the whole domain
# Select one vertical column , t =0 , x=30, y=30
mySounding = sounding() # Create an empty sounding
#Add fields
mySounding.addField('pressure', pressure[0,:,30,30], "Pressure", "Pa")
mySounding.addField('temperature', temperature[0,:,30,30], "Temperature", "C")
mySounding.addField('dewPointTemperature', dewPointTemperature[0,:,30,30], "Dew Point Temperature", "C")
mySounding.quickPlot()
)
主要是利用siphon直接下载数据，未来多种数据包括卫星可能直接绘图!
#https://unidata.github.io/siphon/latest/examples/upperair/Wyoming_Request.html#sphx-glr-download-examples-upperair-wyoming-request-py
from datetime import datetime
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
#from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, SkewT
from metpy.units import units
from siphon.simplewebservice.wyoming import WyomingUpperAir
date = datetime(2017, 9, 10, 6)
station = 'MFL'
df = WyomingUpperAir.request_data(date, station)
print(df.columns)
p = df['pressure'].values * units(df.units['pressure'])
T = df['temperature'].values * units(df.units['temperature'])
Td = df['dewpoint'].values * units(df.units['dewpoint'])
u = df['u_wind'].values * units(df.units['u_wind'])
v = df['v_wind'].values * units(df.units['v_wind'])
fig = plt.figure(figsize=(9, 9))
#add_metpy_logo(fig, 110, 100)
skew = SkewT(fig, rotation=45)
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
# Calculate LCL height and plot as black dot
lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
# Calculate full parcel profile and add to plot as black line
prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)
# Shade areas of CAPE and CIN
#skew.shade_cin(p, T, prof)
#skew.shade_cape(p, T, prof)
# An example of a slanted line at constant T -- in this case the 0
# isotherm
skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
# Show the plot
plt.savefig("asending.png")
plt.show()
# Brian Blaylock
# July 6, 2015
# Download, process, and Plot Sounding Data from University of Wyoming
import urllib3
from bs4 import BeautifulSoup
from skewt import SkewT
from matplotlib import pyplot as plt
stn = ['52203','51777','51839','51828','51709','51644','51463'] #72572 is ID for SLC
name = ['HM','RQ','MF','HT','KS','KQ','Urumqi'] #72572 is ID for SLC
year= '2015'
month = '04'
day = '27'
hour = '12' #either 12 or 00
# 1)
# Wyoming URL to download Sounding from
http = urllib3.PoolManager()
for sta,na in zip(stn,name):
print(sta)
url = 'http://weather.uwyo.edu/cgi-bin/sounding?region=naconf&TYPE=TEXT%3ALIST&YEAR='+year+'&MONTH='+month+'&FROM='+day+hour+'&TO='+day+hour+'&STNM='+sta
print(url)
response = http.request('GET', url)
# 2)
# Remove the html tags
soup = BeautifulSoup(response.data)
data_text = soup.get_text()
#text_file = open(stn+year+month+day+hour+'.txt', "w")
#text_file.write("Purchase Amount: %s" % data_text)
#text_file.close()
#data_text.to_csv(stn+year+month+day+hour+'.csv',encoding='gbk')
# 3)
# Split the content by new line.
splitted = data_text.split("\n",data_text.count("\n"))
# 4)
# Write this splitted text to a .txt document
Sounding_filename = str(sta)+'.'+str(year)+str(month)+str(day)+str(hour)+'.txt'
f = open(Sounding_filename,'w')
for line in splitted[4:]:
f.write(line+'\n')
f.close()
# 5)
print(Sounding_filename)
S = SkewT.Sounding(Sounding_filename)
#sounding.plot_skewt()
#S=SkewT.Sounding("./skewt/examples/bna_day1.txt")
#T=SkewT.Sounding("./skewt/examples/bna_day2.txt")
S.make_skewt_axes()
S.add_profile(color='r',bloc=0)
plt.savefig(sta+year+month+day+hour+'.png')
#S.soundingdata=T.soundingdata # replace the sounding data in S with that from T
#S.add_profile(color='b',bloc=1)
'''
And if you wish to plot more than one sounding on the same plot do this...
S = SkewT.Sounding(filename=Sounding_filename)
T = SkewT.Sounding(filename="72572.2015061212.txt")
S.make_skewt_axes(tmax=55)
S.add_profile(color='r',linewidth=5,bloc=0)
S.soundingdata=T.soundingdata
S. add_profile(color="b",bloc=1)
'''
image.png
IDD
;**************************************************
; skewt_7.ncl
;
; Concepts illustrated:
; - Reading ds336.0 netCDF files
; - Specifying desired station id(s) and extracting desired data
; - Changing non-standard units associated with the time coordinate variable
; - Adding an _FillValue attribute to a variable
; - Looping over files, stations and times
; - Drawing Skew-T plots
; - Creating simple ascii (text) files of each sounding
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
diri = "./" ; directory containing input file(s)
fili = systemfunc("cd "+diri+" ; ls Upperair_20150427*nc") ; all file names
nfili = dimsizes(fili) ; number of files to be used
print(fili)
idkey = (/51839/) ; desired station(s): (/ id1, id2,.../)
nidkey = dimsizes(idkey)
; user names for station(s), one per 'idkey'
idName = (/"Minfeng"/) ; used in plots and text files
txtFile = True ; Is a text file to be created? True=>>Yes, False==>No
txtDir = "./" ; directory for text files if txtFile=True
pltType = "png" ; png, ps, eps, pdf, x11
pltDir = "./" ; directory for plots
;====================> ; end of generic specifications
; skewT options
skewtOpts = True
;skewtOpts@DrawColAreaFill = True ; default is False
dataOpts = True ; options describing data and ploting
dataOpts@HspdHdir = True ; wind speed and dir [else: u,v]
; text file header(s)
varHead = [/" P HGT T TDEW WSPD WDIR "/]
untHead = [/" mb m C C m/s "/]
do nf=0,nfili-1 ; loop over each file
f = addfile(diri+fili(nf), "r")
idnum := f->wmoStaNum
synTime := f->synTime
synTime@units = "seconds since 1970-1-1 00:00:0.0" ; replace: "seconds since (1970-1-1 00:00:0.0)"
ymdh := cd_calendar(synTime, -3)
print(ymdh)
do id=0,nidkey-1 ; loop over each station
ii := ind(idnum.eq.idkey(id)) ; may be multiple times (eg,00,06,12,18)
nii = dimsizes(ii) ; for each station
ii0 = ii(0) ; convenience: use info associated with 1st
staLat = f->staLat(ii0)
staLon = f->staLon(ii0)
staElev = f->staElev(ii0)
staName = tostring(f->staName(ii0,:)) ; char staName(recNum, staNameLen)
do i=0,nii-1 ; loop over time(s) (eg, 00,06,12,18)
; for current station
if (.not.ismissing(ii(i))) then
iii = ii(i) ; convenience
prMan = f->prMan(iii,:) ; 'Man' =" Mandatory levels
wdMan = f->wdMan(iii,:) ; (recNum, manLevel) ==> (:) ==> hPa
wsMan = f->wsMan(iii,:)
htMan = f->htMan(iii,:)
tpMan = f->tpMan(iii,:)
tdMan = f->tdMan(iii,:)
tpMan = tpMan-273.15 ; skewT expects degC
tpMan@units = "degC"
tdMan = tpMan-tdMan
tdMan@units = "degC"
numSigW= f->numSigW(iii) ; 'SigW' = Significant Wind levels
numSigW@_FillValue = 99999 ; _FillValue not on file
if (.not.ismissing(numSigW) .and. numSigW.gt.0) then
dataOpts@PlotWindH = True ; plot wind barbs at height lvls
dataOpts@Height = f->htSigW(iii,:) ; height of wind reports!!!show the height
dataOpts@Hspd = f->wsSigW(iii,:) ; speed [or u component]
dataOpts@Hdir = f->wdSigW(iii,:) ; dir [or v component]
else
dataOpts@PlotWindH = False ; No pibal winds to plot
end if
pltName = "skewt."+idName(id)+"_"+idkey(id)+"."+ymdh(iii) ; eg: skewt.Bermuda_78016.2013030612.png
pltPath = pltDir+pltName
pltType@wkWidth = 800 ; large
pltType@wkHeight = 800
wks = gsn_open_wks (pltType, pltPath)
skewtOpts@tiMainString = idName(id)+": "+idkey(id)+": "+ymdh(iii)
;***********************
;---Draw skewT background
skewtOpts@DrawColLine = True ; draw background lines in color, (False for all black lines)
;skewtOpts@DrawHeightScale = True ; default is False
skewtOpts@DrawFahrenheit = False ; True for deg F, False for deg C
;skewtOpts@tiMainString = "Sounding at "+snd_lat+","+snd_lon+": "+in_file
skewtOpts@tiMainFontHeightF = 0.0155
skewtOpts@vpHeightF = 0.85 ; controls height of plot
skewtOpts@vpWidthF = 0.85 ; controls width of plot
skewtOpts@vpXF = 0.07 ; controls off-set from left
skewtOpts@vpYF = 0.92 ; controls off-set from top
skewtOpts@ThermoInfo =False
;skewtOpts@DrawStandardAtm = False
skewt_bkgd = skewT_BackGround (wks, skewtOpts)
; Draw the skewT plot background with all the overlays
draw (skewt_bkgd)
; Draw the skew-T data on top of everything else
dataOpts = True
dataOpts@colTemperature = "red" ; temperature line color
dataOpts@colDewPt = "blue" ; dewpoint line color
dataOpts@DrawWindBarbThk = 2. ; wind barb thickness
dataOpts@Parcel = 1 ; subscript corresponding to initial parcel
;dataOpts@WspdWdir = False ; True = use wind speed and dir, False = use U & V
dataOpts@Wthin = 0 ; 0 = draw wind barbs at all levels, 1 = draw every other, 2 = every 3rd, etc.
skewt_data = skewT_PlotData (wks, skewt_bkgd, prMan,tpMan,tdMan,htMan \
, wsMan,wdMan, dataOpts)
frame(wks) ; advance frame
print("Created plot: "+pltName+"."+pltType)
if (txtFile) then
varList = [/prMan, htMan, tpMan, tdMan, wsMan, wdMan/]
txtName = idName(id)+"_"+ymdh(iii)+".ManLevels.txt" ; eg: Bermuda_2013030612.ManLevels.txt
txtPath = txtDir+txtName
write_table(txtPath, "w", varHead, "%s") ; "w" => create or overwrite
write_table(txtPath, "a", untHead, "%s") ; "a" => append
write_table(txtPath, "a", varList, "%7.0f%7.0f%7.1f%7.1f%7.1f%7.0f")
end if
end if ; .not.ismissing(ii(i))
end do ; i (time)
end do ; id (station id)
end do ; nf (file)
image.png
模式结果
获取廓线数据按照怀俄明大学变量
from __future__ import print_function, division
from netCDF4 import Dataset
from wrf import getvar, CoordPair, ll_to_xy
import pandas as pd
names=['PRES','HGHT','TEMP','DWPT','RELH','MIXR','DRCT','SKNT','THTA','THTE','THTV']
df = pd.DataFrame(columns = names)
ncfile = Dataset("wrfout_d01_2015-04-27_02:00:00")
x_y = ll_to_xy(ncfile, 41.71,82.95)
print (x_y)
# Get the geopotential height (m) and pressure (hPa).
z = getvar(ncfile, "z")
zp= z[:,x_y[0],x_y[1]]
p = getvar(ncfile, "pressure")
pp= p[:,x_y[0],x_y[1]]
rh = getvar(ncfile, "rh")
rhp= rh[:,x_y[0],x_y[1]]
td = getvar(ncfile, "td")
tdp= td[:,x_y[0],x_y[1]]
tc = getvar(ncfile, "tc")
tcp= tc[:,x_y[0],x_y[1]]
print(tcp)
wspd = getvar(ncfile, "wspd_wdir",units='kt')[0]
wdir = getvar(ncfile, "wspd_wdir",units='kt')[1]
wspdp= wspd[:,x_y[0],x_y[1]]
wdirp= wdir[:,x_y[0],x_y[1]]
df.PRES=pp
df.HGHT=zp
df.RELH=rhp
df.DWPT=tdp
df.TEMP=tcp
df.DRCT=wdirp
df.SKNT=wspdp
df.MIXR=None
df.THTA=None
df.THTE=None
df.THTV=None
df.to_csv('Resulta.txt',sep=' ',float_format='%.2f',header=True, index=False,na_rep='NaN')
;**************************************************
; skewt_10.ncl
;
; Concepts illustrated:
; - Read from WRF netCDF and extracting necessary variables
; - Using 'wrf_user_ll_to_ij' get indices (subscripts) near specified location
; - Using NCL's 'wind_component' to convert wind speed and direvtion to u,v
; - Adjusting the fortran 1-based subscripts to NCL 0-based subscripts
; - Plotting a hodograph onto a skew-T plot
;**************************************************
; Author: Joesph Grim
; Project Scientist
; Aviation Applications Program
; Research Application Laboratory
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
begin
;---Set parameters
snd_lat = 39.33 ; latitude of model profile
snd_lon = 89.42 ; longitude of model profile
hodo_top_pres = 150. ; top of hodograph in hPa
in_dir = "/public/home/hysplit/work/WRF381/WRF-381/test/em_real/"
in_file = "wrfout_d01_2015-04-30_00:00:00" ; WRF output file
;---Read in WRF data
in_path = in_dir + in_file
ncdf_in = addfile(in_path+".nc","r") ; add .nc extension
P_tot = wrf_user_getvar(ncdf_in,"pressure",0) ; model 3D pressure
z_tot = wrf_user_getvar(ncdf_in,"z",-1) ; model 3D height
uvmet = wrf_user_getvar(ncdf_in,"uvmet",-1) ; Earth-relative U and V winds on mass points
printVarSummary(uvmet)
U = uvmet(0,0,:,:,:) ; Earth-relative U wind on mass points
V = uvmet(1,0,:,:,:) ; Earth-relative V wind on mass points
;wspd = wind_speed(U,V)
wspd =sqrt(U^2+V^2)
wdir = wind_direction(U,V,0)
TC = wrf_user_getvar(ncdf_in,"tc",-1) ; model temperature in Celsius
TD = wrf_user_getvar(ncdf_in,"td",-1) ; model dewpoint in Celsius
;---Determine the model i,j indices for the given lat/lon
loc = wrf_user_ll_to_ij(ncdf_in,snd_lon,snd_lat,True) ; model i,j indices for the given lat/lon
locX = loc(0) - 1 ; subtract 1, since WRF is base 1, but NCL is base 0
locY = loc(1) - 1 ; subtract 1, since WRF is base 1, but NCL is base 0
;***********************
;---Begin to create plot
;***********************
PlotType = "png"
;PlotName = "skewT_hodograph"
PlotName = "skewt"
PlotType@wkWidth = 800 ; large
PlotType@wkHeight = 800
wks = gsn_open_wks(PlotType,PlotName) ; open workstation
;---Draw skewT background
skewtOpts = True
skewtOpts@DrawColLine = True ; draw background lines in color, (False for all black lines)
skewtOpts@DrawFahrenheit = False ; True for deg F, False for deg C
skewtOpts@tiMainString = "Sounding at "+snd_lat+","+snd_lon+": "+in_file
skewtOpts@tiMainFontHeightF = 0.0155
skewtOpts@vpHeightF = 0.85 ; controls height of plot
skewtOpts@vpWidthF = 0.85 ; controls width of plot
skewtOpts@vpXF = 0.07 ; controls off-set from left
skewtOpts@vpYF = 0.92 ; controls off-set from top
skewt_bkgd = skewT_BackGround (wks, skewtOpts)
; Draw the skewT plot background with all the overlays
draw (skewt_bkgd)
; Draw the skew-T data on top of everything else
dataOpts = True
dataOpts@PlotWindH = True ; plot wind barbs at height lvls
dataOpts@Height = z_tot ; height of wind reports
dataOpts@Hspd = wspd ; speed [or u component]
dataOpts@Hdir = wdir ; dir [or v component]
dataOpts@colTemperature = "red" ; temperature line color
dataOpts@colDewPt = "blue" ; dewpoint line color
dataOpts@DrawWindBarbThk = 2. ; wind barb thickness
dataOpts@Parcel = 1 ; subscript corresponding to initial parcel
dataOpts@WspdWdir = True ; True = use wind speed and dir, False = use U & V
dataOpts@Wthin = 0 ; 0 = draw wind barbs at all levels, 1 = draw every other, 2 = every 3rd, etc.
skewT_data = skewT_PlotData(wks, skewt_bkgd, P_tot(:,locY,locX), \
TC(0,:,locY,locX), \
TD(0,:,locY,locX), \
z_tot(0,:,locY,locX), \
wspd(:,locY,locX), \
wdir(:,locY,locX), \
dataOpts)
frame(wks) ; advance frame
print("Created plot: "+PlotName+"."+PlotType)
end
读取python输出的txt
;**************************************************
; skewt_2.ncl
;
; Concepts illustrated:
; - Drawing Skew-T plots
; - Drawing a Skew-T plot with and without wind barbs
; - Customizing the background of a Skew-T plot
;*************************************************
;
; This file is loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;
; This file still has to be loaded manually
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
;***********************************************
begin
; --- Read Data; no winds; set to _FillValue------------;
diri = "./"
fili = "sounding.testdata22"
nlvl = 36
ncol = 11
TestData = asciiread (diri+fili , (/nlvl,ncol/), "float")
p = TestData (:,0)
z = TestData (:,1)
;q = TestData (:,4)
tc = TestData (:,2) + 2. ; for demo purposes
tdc = TestData (:,3)
wspd = TestData (:,7)
wdir = TestData (:,6)
; create BOGUS winds
;wspd = new ( nlvl, "float", -999.) ; wspd@_FillValue = -999.
;wdir = new ( nlvl, "float", -999.) ; wdir@_FillValue = -999.
;*************************
; create plot
;*************************
wks = gsn_open_wks ("png", "skewt") ; send graphics to PNG file
; --- Create background skew-T; plot sounding ---------------------
skewtOpts = True
skewtOpts@DrawColAreaFill = True ; default is False
skewtOpts@tiMainString = "Raob Data; No Winds"
dataOpts = False ; no options
; dataOpts@PrintZ = False
skewt_bkgd = skewT_BackGround (wks, skewtOpts)
skewt_data = skewT_PlotData (wks, skewt_bkgd, p,tc,tdc,z \
, wspd,wdir, dataOpts)
draw (skewt_bkgd)
draw (skewt_data)
frame(wks)
; --- Create background skew-T and plot sounding + winds----------------
; Create winds for demo
wspd = fspan (0., 150., nlvl) ; wind speed at each level
wdir = fspan (0., 360., nlvl) ; wind direction
; Create a few "pibal" reports
hght = (/ 1000., 3000., 7000., 25000. /)/3.208 ; hgt in M
hspd = (/ 50., 27., 123., 13. /) ;speed at each height
hdir = (/ 95., 185., 275., 355. /) ;direction
dataOpts = True ; options describing data and ploting
dataOpts@PlotWindH = True ; if available, plot wind barbs at height lvls
dataOpts@HspdHdir = True ; wind speed and dir [else: u,v]
dataOpts@Height = hght ; height of wind reports
dataOpts@Hspd = hspd ; speed [or u component]
dataOpts@Hdir = hdir ; dir [or v component]
skewtOpts@tiMainString = "Raob; [Wind Reports]"
skewtOpts@DrawColAreaFill = True ; default is False
skewtOpts = True
skewtOpts@DrawColAreaFill = True ; default is False
skewt_bkgd = skewT_BackGround (wks, skewtOpts)
skewt_data = skewT_PlotData (wks, skewt_bkgd, p,tc,tdc,z \
, wspd,wdir, dataOpts)
draw (skewt_bkgd)
draw (skewt_data)
frame(wks)
; --- Make the wind barbs thicker (NCL 6.4.0)
skewtOpts@tiMainString = "Raob; Thicker Wind Barbs"
dataOpts@DrawWindBarbThk = 1.5 ; default value is 1.0 (may be fractional)
skewt_bkgd = skewT_BackGround (wks, skewtOpts)
skewt_data = skewT_PlotData (wks, skewt_bkgd, p,tc,tdc,z \
, wspd,wdir, dataOpts)
draw (skewt_bkgd)
draw (skewt_data)
frame(wks)
end
image.png
