Hi,
1)triyng to get temperature data from gfs opendap server but script download the from a wrong latitude / longitude
the correct one are 45N and 9E (north italy)
Ugo
script
----------------------------------------------------------------------
; This example shows how to read geographic data
; from Natural Earth shapefiles
; and plot them as polylines and polygons.
;----------------------------------------------------------------------
; This particular example plots data for Switzerland.
;----------------------------------------------------------------------
; Download the shapefiles from http://www.naturalearthdata.com/
; Unzip to a directory
;----------------------------------------------------------------------
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/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
;******************
; read gfs data *
;******************
url = "http://nomad1.ncep.noaa.gov:9090/dods/gfs_master/gfs" + systemfunc("date +%Y%m%d") + "/"
systemdate = systemfunc("date +%H")
if ((systemdate.ge.06).and.(systemdate.lt.12)) then
filename = url + "gfs_master_00z"
end if
if ((systemdate.ge.12).and.(systemdate.lt.18)) then
filename = url + "gfs_master_00z"
end if
if ((systemdate.ge.18).and.(systemdate.lt.24)) then
filename = url + "gfs_master_06z"
end if
if ((systemdate.ge.24).and.(systemdate.lt.06)) then
filename = url + "gfs_master_18z"
end if
exists = isfilepresent(filename)
if(.not.exists) then
print("OPeNDAP isfilepresent test unsuccessful.")
print("Either file doesn't exist, or NCL does not have OPeNDAP capabilities on this system")
else
print("OPeNDAP isfilepresent test successful.")
gfs = addfile(filename,"r")
vars = getfilevarnames(gfs)
end if
;---data for 2 metres temperature plot
NTIMES = dimsizes(gfs->time) ; number of times in the file
t=new((/NTIMES,2/), string)
s=new((/1,NTIMES/), string)
do it1 = 0,NTIMES-1
;do it2 = 0,1
t(it1,0)= it1
print(t(it1,0))
;print(TEMP2M)
t(it1,1) = (gfs->tmp2m(it1,45,9)) - 273.15
print(t(it1,1))
s(0,it1) = t(it1,0)+" "+t(it1,1)
print(gfs->tmp2m(it1,45,9))
;end do
end do
;end do
s=(/ s(0,ispan(0,NTIMES-1,1)) /)
asciiwrite("prova.txt", s)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sat Sep 17 11:31:32 2011
This archive was generated by hypermail 2.1.8 : Thu Sep 22 2011 - 17:12:44 MDT