[ncl-talk] Problem reading data from NCEP DODS server
Mary Haley
haley at ucar.edu
Thu May 19 10:42:41 MDT 2016
Hi Chris,
Whenever you do a calculation with NCL and assign it to a new variable,
this causes all metadata, except for the _FillValue attribute (if any), to
be stripped off.
For example, if "tempC" is a variable containing temperatures in degrees
celsius, and it has a "long_name" attribute of "temperature" and a "units"
of "degC", then doing this:
tempF = 1.8 * tempC + 32.
is going to cause "tempF" to have no attributes at all. When you do this,
however:
tempF = tempC
tempF = 1.8 * tempC + 32.
then "tempF" is first an exact copy of "tempC", and hence will have all of
its metadata attached. When you do the calculation, the metadata will not
be stripped from the existing "tempF".
As you pointed out, this is not ideal because now the "units" and
potentially other attributes are wrong.
But, if this is the only thing that's wrong, then you can easily fix it
with:
tempF = tempC
tempF = 1.8 * tempC + 32.
tempF at units = "degF" ; Fix the units
This is perfectly okay in NCL, but I agree that it's a little odd.
The other way around this is to copy the metadata yourself using
"copy_VarMeta":
tempF = 1.8 * tempC + 32. ; metadata is stripped
copy_VarMeta(tempC,tempF) ; copies *all* metadata from tempC to tempF
tempF at units = "degF" ; fix the units
If you just want to copy the variable coordinate arrays, but not the
variable attributes, then:
tempF = 1.8 * tempC + 32.
copy_VarCoords(tempC,tempF) ; copies only the coordinate arrays (if any)
from tempC to tempF
;---You wil have to add all the attributes yourself
tempF at long_name = "temperature
tempF at units = "degF" ; still need to fix the units
See:
http://www.ncl.ucar.edu/Document/Functions/Contributed/copy_VarMeta.shtml
http://www.ncl.ucar.edu/Document/Functions/Contributed/copy_VarCoords.shtml
http://www.ncl.ucar.edu/Document/Functions/Contributed/copy_VarAtts.shtml
--Mary
On Wed, May 18, 2016 at 7:09 PM, Herbster, Christopher G. <herbstec at erau.edu
> wrote:
> 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","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> 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 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> 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 Office
>>
>> 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
>> 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/70d38060/attachment.html
More information about the ncl-talk
mailing list