[ncl-talk] Problem reading data from NCEP DODS server
Herbster, Christopher G.
herbstec at erau.edu
Wed May 18 19:09:13 MDT 2016
Mary,
Thanks again for your help. First, I forgot to mention that the nomads server only holds two days of data. The date in the script had the value hardwired to a date that is no longer there. Just change it to "yesterday" and it should find data on the other end.
Next, I was getting error/warnings about not having coordinates for a variable I was plotting. When I sent my message, I didn't really know where that was coming from. I now know that this was from the "speed" variable that I created from the U and V components of the current. I got rid of this issue by what I would call a hammer technique. I first did
speed = u_velocity
and then
speed = sqrt((u_velocity^2)+(v_velocity^2)) * 1.9438 ; convert m/s to knots
This got rid of the error, and, I think, has also fixed any plotting issues there may have been with presenting the data on a map.
However, this hammer is more than I wanted, as now all of the attributes are set to the U component of the current. Is there a better way to do this? Should I do this and then change just what I need to change, or should I just create the code to set all the metadata (attributes?) for the variable?
I'd like to start adopting a "best practice" method, as I hope to eventually have my students develop code in NCL for some of my class projects.
Thanks again for your help, and sorry for the cryptic questions. I'll try to do better. (Though I think I've gotten some insight on how to deal with my student's questions! (-:)
Cheers,
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 <haley at ucar.edu>
Sent: Wednesday, May 18, 2016 11:06 AM
To: Herbster, Christopher G.
Cc: ncl-talk at ucar.edu
Subject: Re: [ncl-talk] Problem reading data from NCEP DODS server
Chris,
Yes, the "all(ismissing(sst))" is a way to trap this special case.
I'm not sure what you are asking in your second question. When you subset "sst" with the given lat/lon range, then the lat/lon coordinate arrays should also be subsetted automatically.
I just tried to run your script again to verify that you indeed have coordinate arrays attached to your data, but it appears the file is no longer there.
You can verify this yourself by doing a "printVarSummary" on "sst" after you read it off a file.
Here's a little sample script that you should be able to run as-is:
a = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/uv300.nc<http://uv300.nc>","r")
u = a->U(1,:,:)
printVarSummary(u)
u_subset = a->U(1,{0:30},{-10:10})
printVarSummary(u_subset)
The output (trimmed) will look like this:
Variable: u
Type: float
. . .
Dimensions and sizes: [lat | 64] x [lon | 128]
Coordinates:
lat: [-87.8638..87.8638]
lon: [-180..177.1875]
Number Of Attributes: 5
time : 7
_FillValue : -999
long_name : Zonal Wind
short_name : U
units : m/s
Variable: u_subset
Type: float
. . .
Dimensions and sizes: [lat | 11] x [lon | 7]
Coordinates:
lat: [1.395307..29.30136]
lon: [-8.4375..8.4375]
Number Of Attributes: 5
time : 7
_FillValue : -999
long_name : Zonal Wind
short_name : U
units : m/s
Note that for "u_subset", the lat/lon coordinates have the new subsetted range, based on the closest *less than* values to the lat/lon range I gave it.
--Mary
On Tue, May 17, 2016 at 5:54 PM, Herbster, Christopher G. <herbstec at erau.edu<mailto:herbstec at erau.edu>> wrote:
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<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/
From: Mary Haley [mailto:haley at ucar.edu<mailto:haley at ucar.edu>]
Sent: Tuesday, May 17, 2016 5:58 PM
To: Herbster, Christopher G. <herbstec at erau.edu<mailto:herbstec at erau.edu>>
Cc: ncl-talk at ucar.edu<mailto: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/20160519/e8abc17a/attachment.html
More information about the ncl-talk
mailing list