Gfs opendap wrong latitude / longitude data

From: ugo merlini <ugomerlini_at_nyahnyahspammersnyahnyah>
Date: Sat Sep 17 2011 - 09:31:24 MDT

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