[ncl-talk] Problem reading data from NCEP DODS server

Herbster, Christopher G. herbstec at erau.edu
Tue May 17 17:54:36 MDT 2016


Mary,

Thanks for checking what seems obvious in hindsight.  Feeling a bit sheepish.

Can I trap this with an:

if(all(ismissing(sst))) then
 … ; don’t do much
else
… ; make plots
end if

Boy have I wasted some time thinking this had to do with the longitude being 74 – 434 degrees east.  I had originally thought that I was somehow only getting data over land, where SST would certainly be missing.

Also, is there a preferred way to have the meta data for the subset data have ranges that are indicative of the subset, rather than the whole domain?

Thanks again for catching my oversight,

Chris H.
--

Dr. Christopher G. Herbster
Associate Professor
Director of Science and Technology
for the ERAU Weather Center
Applied Aviation Sciences
Embry-Riddle Aeronautical Univ.
600 S. Clyde Morris Blvd.
Daytona Beach, FL 32114-3900
 386.226.6444 Office
386.226.6446 Weather Center
http://wx.erau.edu/

Schedule at:  http://wx.erau.edu/faculty/herbster/Schedules/

From: Mary Haley [mailto:haley at ucar.edu]
Sent: Tuesday, May 17, 2016 5:58 PM
To: Herbster, Christopher G. <herbstec at erau.edu>
Cc: ncl-talk at ucar.edu
Subject: Re: [ncl-talk] Problem reading data from NCEP DODS server

Chris,

I think the issue is simply that the first time step of "sst" is all missing.

If you select a different time step, then it looks fine:

; sst = f->sst(mytime,:,Slat:Nlat,Wlon:Elon)

sst = f->sst(:,:,Slat:Nlat,Wlon:Elon)    ; NOTE: READING ALL TIME STEPS
print("done reading sst")
print("dimsizes sst = " + str_join(dimsizes(sst),","))

do nt=0,dimsizes(sst(:,0,0,0))-1
  printMinMax(sst(nt,:,:,:),0)
end do

This is the output:

(0)  dimsizes sst = 65,1,265,481
nt = 0
sea_surface_temperature (c)  : min=1.26765e+30   max=1.26765e+30
nt = 1
sea_surface_temperature (c)  : min=21.6418   max=32.2542
nt = 2
sea_surface_temperature (c)  : min=21.6017   max=32.3401
nt = 3
sea_surface_temperature (c)  : min=21.5984   max=32.3429
nt = 4
sea_surface_temperature (c)  : min=21.6628   max=32.1836
nt = 5
sea_surface_temperature (c)  : min=22.0586   max=32.491
nt = 6
sea_surface_temperature (c)  : min=22.2371   max=32.6994

--Mary


On Tue, May 17, 2016 at 10:28 AM, Herbster, Christopher G. <herbstec at erau.edu<mailto:herbstec at erau.edu>> wrote:
Hi folks,

I’m trying to retrieve some data for a subset region from the Global Real Time Ocean Forecast System (RTOFS), and I can’t figure out why I get all missing values when I try to read the data.

This seems to be the case for me even if I try to read this whole domain.

Here is some code:
------
; read-rtofs-dods.ncl
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
;

Dev = "png"

;
myDate = 20160516

;NOMADS = "http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global"+myDate+"/rtofs_glo_2ds_forecast_daily_prog"

NOMADS = "http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global"+myDate+"/rtofs_glo_2ds_forecast_3hrly_prog"

filename = NOMADS ; = url
print("Fetching URL: ")
print(" "+filename)

  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")
    exit
  else
    print("OPeNDAP isfilepresent test successful.")
    f = addfile(filename,"r")
    vars = getfilevarnames(f)
    print(vars)   ; should be (in any order):

    if(.not.any(ismissing(vars))) then
      do i=0,dimsizes(vars)-1
        printFileVarSummary (f,vars(i))
      end do
    end if
  end if

; Define boundary of data
; Assign lat/lon boundaries to the data
; Longitudes between 0 and 74 E need to have 360 added to them
;
myNlat = 32.0
mySlat = 10.0
myWlon = 260.0
myElon = 300.0

;myNlat = 26.0
;mySlat = 21.0
;myWlon = 260.0
;myElon = 300.0

; Full domain
;myNlat = 90.0
;mySlat = -90.0
;myWlon = 74.0
;myElon = 434.0



; Get lat/lon/time info first
lat = f->lat
print ("Done reading lat")
lon = f->lon
print ("Done reading lon")
time = f->time(0)
print ("Done reading time")
;
; Use just one time when "0"
mytime = 0
;mytime = ":"

Nlat = closest_val(myNlat, lat)
print ("Nlat = " + Nlat + " Lat = " + lat(Nlat) )

Slat = closest_val(mySlat, lat)
print ("Slat = " + Slat + " Lat = " + lat(Slat) )

Wlon = closest_val(myWlon, lon)
print("Wlon = " + Wlon + " Lon = " + lon(Wlon) )

Elon = closest_val(myElon, lon)
print("Elon = " + Elon + " Lon = " + lon(Elon) )
;
;exit
;
mylat = lat(Slat:Nlat)
mylon = lon(Wlon:Elon)

printVarSummary (mylat)
printVarSummary (mylon)

print (mylat)

;exit


delete(lat)
delete(lon)

lat = f->lat(Slat:Nlat)
lon = f->lon(Wlon:Elon)
print ("Done reading subset lat/lon")

printVarSummary(lat)
print (lat)

printVarSummary(lon)
print (lon)

;exit


sst = f->sst(mytime,:,Slat:Nlat,Wlon:Elon)
;sst = f->sst
print ("Done reading sst")
printVarSummary(sst)
print (sst)

; exit  ; Uncomment to stop after sst data are read

u_velocity = f->u_velocity(mytime,:,Slat:Nlat,Wlon:Elon)
;u_velocity = f->u_velocity
print ("Done reading U")
printVarSummary(u_velocity)
print (u_velocity)
v_velocity = f->v_velocity(mytime,:,Slat:Nlat,Wlon:Elon)
;v_velocity = f->v_velocity
print ("Done reading V")
printVarSummary(v_velocity)
print (v_velocity)
speed = sqrt((u_velocity*u_velocity) + (v_velocity*v_velocity))
print ("Done calculating Speed")
print (speed)


exit
;;;;;;
;;;;;;  END of code snippet

------

If you run this code, I suggest a redirect into a file, as there are a lot of lines written to standard out.

What I don’t understand is that I can read a subset of the lat and lon data, but when I try to get any other data they all come up as missing/fill values.

Any help is GREATLY appreciated and I have really chased my tail trying to get past this.

Thanks,
Chris Herbster
--

Dr. Christopher G. Herbster
Associate Professor
Director of Science and Technology
for the ERAU Weather Center
Applied Aviation Sciences
Embry-Riddle Aeronautical Univ.
600 S. Clyde Morris Blvd.
Daytona Beach, FL 32114-3900
 386.226.6444<tel:386.226.6444> Office
386.226.6446<tel:386.226.6446> Weather Center
http://wx.erau.edu/

Schedule at:  http://wx.erau.edu/faculty/herbster/Schedules/


_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu<mailto:ncl-talk at ucar.edu>
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160517/9e47a936/attachment-0001.html 


More information about the ncl-talk mailing list