Re: Wind speed problem with no wind

Division by 0.0 ... well, then it is up to *you* to do something.

  ugrd = gfs->ugrd_4572m
  printVarSummary(ugrd)

  uu = where(ugrd.eq.0.0, ugrd@_FillValue, ugrd)
  VELOCITA_VENTO =
  round(((gfs->vgrd_4572m)/sin(atan2((gfs->vgrd_4572m),uu)))*3.6,3)

On 10/27/2011 01:25 PM, ugo merlini wrote:
> Hi,
>
> the problem in that U V variables can't assume 0 as value because the line
>
> VELOCITA_VENTO (wind speed) =
> round(((gfs->vgrd_4572m)/sin(atan2((gfs->vgrd_4572m),(gfs->ugrd_4572m))))*3.6,3)
>
> give me a fatal error (div by 0 error)
>
> and missing value for these variables is 9.999E20
>
> I try with where function with no results
>
> NLATITUDE= dimsizes(gfs->lat)
> NLONGITUDE= dimsizes(gfs->lon)
> do latitudine = 0,NLATITUDE-1
> do longitudine = 0,NLONGITUDE-1
> vento_velocita_plot =
> gsn_csm_contour(wks,where((gfs->ugrd_4572m(periodo,latitudine,longitudine).gt.0).and.(gfs->vgrd_4572m(periodo,latitudine,longitudine).gt.0).and.(.not.ismissing(gfs->ugrd_4572m(periodo,latitudine,longitudine))).and.(.not.ismissing(gfs->vgrd_4572m(periodo,latitudine,longitudine)))
> ,
> round(((gfs->vgrd_4572m(periodo,latitudine,longitudine))/sin(atan2((gfs->vgrd_4572m(periodo,latitudine,longitudine)),(gfs->ugrd_4572m(periodo,latitudine,longitudine)))))*3.6,3),0),vento_velocita)
> end do
> end do
> Regards
>
> Ugo
>
> ------------------------------------------------------------------------
> Subject: Re: Wind speed problem with no wind
> From: haley@ucar.edu
> Date: Mon, 24 Oct 2011 09:13:38 -0600
> CC: ncl-talk@ucar.edu
> To: ugomerlini@hotmail.com
>
> Ugo,
>
> It is too difficult to debug such a big piece of code with out the data
> files to run it.
>
> I think the problem might be with this code:
>
> U = gfs->ugrd_4572m * 3.6
> V = gfs->vgrd_4572m * 3.6 ; Convective rain
>
> When you do a calculation like this, all metadata is stripped from the
> file. To prevent this, try:
>
> ;---Metadata is retained here
> U = gfs->ugrd_4572m
> V = gfs->vgrd_4572m
>
> ;---Now convert to different units
> U = U * 3.6
> V = V * 3.6
>
> If U and V contain a "units" attribute, you should update it:
>
> U@units = "..." ; replace ... with correct units
> V@units = "..."
>
> You can debug your code yourself using printVarSummary to help determine
> if your variables contain the necessary coordinate information:
>
> printVarSummary(U)
> printVarSummary(V)
>
> The above tells you a lot about the data: the type, size, whether it has
> coordinate arrays, what attributes it contains, and whether it has a
> missing value ("_FillValue" attribute).
>
> In the future, please try debugging the code yourself, or provide us
> with data files.
>
> --Mary
>
> On Oct 23, 2011, at 11:00 AM, ugo merlini wrote:
>
> Hi,
>
> I still get some problem. This time is the plotting of wind speed
> where data contain Uwind or Vwind is 0 (no wind) as value. I'm able
> to pass this problem for a specific point but not for an area.
>
> can someone give me some Hints?
>
> Regards
>
> Ugo
>
>
>
> 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"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
>
>
> ;****************************
> ; Funzioni personalizzate *
> ;****************************
>
> ; funzione Data validità
>
> undef("Nuova_data")
> function Nuova_data(data:string)
> local anno, mese, giorno, ora, formato_data, data,
> data_intermedia_1, data_intermedia_2, data_intermedia_3,
> data_intermedia_4, data_intermedia_5, data_intermedia_6
>
> begin
>
> anno = stringtointeger( str_get_cols( data, 0, 3) )
> mese = stringtointeger( str_get_cols( data, 4, 5) )
> giorno = stringtointeger( str_get_cols( data, 6, 7) )
> ora = stringtointeger( str_get_cols( data, 8, 9) )
> formato_data = "hours since 1800-01-01 00:00:00"
>
> data_intermedia_1 = cd_inv_calendar( anno, mese, giorno, ora, 0, 0,
> formato_data, 0 )
>
> ; since data_intermedia_1 has the units "hours since ..." the 2 here
> adds 2 hours
> data_intermedia_2 = data_intermedia_1 + 2
> data_intermedia_2@units = formato_data
>
>
> data_intermedia_3 = cd_calendar( data_intermedia_2, -3 )
>
>
> data_intermedia_4 = sprinti ( "%0.8i", data_intermedia_3 )
>
> ymd2String=str_get_cols ( data_intermedia_4, 0, 7 )
> hr2String =str_get_cols ( data_intermedia_4, 8, 9 )
>
> data_intermedia_5 = systemfunc( "date -d '"+ymd2String+"
> "+hr2String+"' '+%A, %d %B %Y ore %H:00'" )
>
>
>
>
> return (data_intermedia_5)
>
> end
>
>
> ;**********************
> ; Codice principale *
> ;**********************
>
> begin
>
>
> ;*****************************************
> ; Open workstation and define colormap *
> ;*****************************************
>
> wks_type = "png" ; use "newpdf" instead of "pdf" for smaller files
> wks_type@wkWidth = 1024
> wks_type@wkHeight = 1024
> wks = gsn_open_wks(wks_type,"Vento_file_temporaneo")
> gsn_define_colormap(wks,"rainbow+white+gray") ; choose colormap
>
>
> ;******************
> ; read gfs data *
> ;******************
>
> systemdate = systemfunc("date +%H%M")
>
> if ((systemdate.eq."0645"))then
> url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" +
> systemfunc("date +%Y%m%d") + "/"
> filename = url + "gfs_hd_00z"
> run = "Run GFS 00 del " + systemfunc("date +%d/%m/%Y")
> run_nome_file = "Run_GFS_00Z"
> end if
> if ((systemdate.eq."1245")) then
> url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" +
> systemfunc("date +%Y%m%d") + "/"
> filename = url + "gfs_hd_06z"
> run = "Run GFS 06 del " + systemfunc("date +%d/%m/%Y")
> run_nome_file = "Run_GFS_06Z"
> end if
> if ((systemdate.eq."1845")) then
> url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" +
> systemfunc("date +%Y%m%d") + "/"
> filename = url + "gfs_hd_12z"
> run = "Run GFS 12 del " + systemfunc("date +%d/%m/%Y")
> run_nome_file = "Run_GFS_12Z"
> end if
> if ((systemdate.eq."0045")) then
> url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" +
> systemfunc("date -d '-1 day' +'%Y%m%d'") + "/"
> filename = url + "gfs_hd_18z"
> run = "Run GFS 18 del " + systemfunc("date -d '-1 day' +'%d/%m/&Y'")
> run_nome_file = "Run_GFS_18Z"
> end if
>
> url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20111022/"
> filename = url + "gfs_hd_06z"
> run = "Run GFS 06 del " + systemfunc("date +%d/%m/%Y")
> run_nome_file = "Run_GFS_00Z"
>
> 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
>
>
> ;*********************************
> ; risorse per la mappa di base *
> ;*********************************
>
> mappa = True
> mappa@mpProjection = "LambertConformal" ; proiezione mappa
>
> mappa@gsnMaximize = True ; massimizza la grafica in base alle
> dimensioni del foglio
> mappa@gsnDraw = False ; la grafica non viene disegnata in automatico
> mappa@gsnFrame = False ; il frame non avanza in automatico
>
> mappa@mpLimitMode = "LatLon" ; tipo di limite della mappa
> mappa@mpMinLonF = -43.0 ; estremo est della mappa
> mappa@mpMaxLonF = 47.0 ; estremo ovest della mappa
> mappa@mpMaxLatF = 88.5 ; estremo nord della mappa
> mappa@mpMinLatF = 27.0 ; estremo sud della mappa
>
> mappa@mpGridAndLimbOn = True ; disegna i meridiani e i paralleli
> sulla mappa
> mappa@mpGridLineDashPattern = 10 ; tipo di linea utilizzata per
> disegnare i meridiani e i paralleli sulla mappa
>
>
> ;******************************
> ; risorse generali dei plot *
> ;******************************
>
> plot_generale= True
> plot_generale@gsnMaximize = True ; massimizza la grafica in base
> alle dimensioni del foglio
> plot_generale@gsnDraw = False ; la grafica non viene disegnata in
> automatico
> plot_generale@gsnFrame = False ; il frame non avanza in automatico
> plot_generale@gsnAddCyclic = True
> plot_generale@gsnSpreadColors = False ; prendo possesso della
> colormap in modo da personalizzare i colori con i valori delle
> singole mappe
> plot_generale@gsnRightString = "" ; elimino il testo standard
> scritto in alto a destra
> plot_generale@gsnLeftString = "" ; elimino il testo standard scritto
> in alto a sinistra
>
>
> ;*******************************
> ; risorse generali dei testi *
> ;*******************************
>
> testo_generale = True
> testo_generale@txFont = "times-roman" ; stalisco il carattere da
> utilizzare nei testi
> testo_generale@txBackgroundFillColor = "Transparent" ; stalisco il
> colore di sfondo del testo
>
> testo_titolo = testo_generale
> testo_titolo@txFontHeightF = 0.02 ; stabilisco la grandezza del
> testo relativo al titolo
> testo_titolo@txJust = "CenterLeft" ; stabilisco la giustificazione
> del testo relativo al titolo
> testo_titolo@txPosXF = 0.02 ; stabilisco la partenza sulla griglia
> orizontale del testo relativo al titolo
> testo_titolo@txPosYF = 0.885 ; stabilisco la partenza sulla griglia
> verticale del testo relativo al titolo
>
> testo_modello = testo_generale
> testo_modello@txFontHeightF = 0.015 ; stabilisco la grandezza del
> testo relativo al nome del modello
> testo_modello@txJust = "BottomLeft" ; stabilisco la giustificazione
> del testo relativo al nome del modello
> testo_modello@txPosXF = 0.02 ; stabilisco la partenza sulla griglia
> orizontale del testo relativo al nome del modello
> testo_modello@txPosYF = 0.851 ; stabilisco la partenza sulla griglia
> verticale del testo relativo al nome del modello
>
> testo_data = testo_generale
> testo_data@txFontHeightF = 0.015 ; stabilisco la grandezza del testo
> relativo al periodo di riferimento
> testo_data@txJust = "BottomRight" ; stabilisco la giustificazione
> del testo relativo al periodo di riferimento
> testo_data@txPosXF = 0.98 ; stabilisco la partenza sulla griglia
> orizontale del testo relativo al periodo di riferimento
> testo_data@txPosYF = 0.85 ; stabilisco la partenza sulla griglia
> verticale del testo relativo al periodo di riferimento
>
>
> ;******************************
> ; generazione mappa di base *
> ;******************************
>
> map = gsn_csm_map(wks,mappa)
>
> fnames = "/mnt/internetserver/map/shapefile/europa/aaa_full/" +
> (/"europa"/) + ".shp" ; Files to be open
> linecolors = (/"Black"/) ; color definition for each shapefile
> fillcolors = (/"Transparent"/) ; Fill color definition for each
> shapefile
> thicks = (/2/) ; Thickness for each shapefile
> lnres = True ; resources for polylines
> plres = True
> prims = True
> lines = True
> do n=0,dimsizes(fnames)-1 ; Loop through files that we want to read
> geographic information from.
> f = addfile(fnames(n),"r") ; Open the shapefile number n.
> segments = f->segments ; Read data off the shapefile
> geometry = f->geometry
> segsDims = dimsizes(segments)
> geomDims = dimsizes(geometry)
> geom_segIndex = f@geom_segIndex ; Read global attributes
> geom_numSegs = f@geom_numSegs
> segs_xyzIndex = f@segs_xyzIndex
> segs_numPnts = f@segs_numPnts
> geometry_type = f@geometry_type
> numFeatures = geomDims(0)
> lon = f->x
> lat = f->y
> if (geometry_type.eq."polygon") then ; Put if statement outside the loop
> plres@gsFillColor = fillcolors(n)
> plres@gsEdgesOn = True ; Draw border around polygons
> plres@gsEdgeColor = linecolors(n)
> plres@gsEdgeThicknessF = thicks(n)
> do i=0, numFeatures-1 ; Section to draw polygons on map.
> startSegment = geometry(i, geom_segIndex)
> numSegments = geometry(i, geom_numSegs)
> do seg=startSegment, startSegment+numSegments-1
> startPT = segments(seg, segs_xyzIndex)
> endPT = startPT + segments(seg, segs_numPnts) - 1
> dumstr = unique_string("lines") ; This call adds the polygon.
> map@$dumstr$ = gsn_add_polygon(wks, map , lon(startPT:endPT),
> lat(startPT:endPT), plres)
> end do
> end do
> else
> lnres@gsLineThicknessF = thicks(n)
> lnres@gsLineColor = linecolors(n)
> do i=0, numFeatures-1 ; Section to draw polylines on map.
> startSegment = geometry(i, geom_segIndex)
> numSegments = geometry(i, geom_numSegs)
> do seg=startSegment, startSegment+numSegments-1
> startPT = segments(seg, segs_xyzIndex)
> endPT = startPT + segments(seg, segs_numPnts) - 1
> dumstr = unique_string("primitive") ; This call adds the line segment.
> map@$dumstr$ = gsn_add_polyline(wks, map, lon(startPT:endPT),
> lat(startPT:endPT), lnres)
> end do
> end do
> end if
> delete(lat) ; Clean up before we read in same variables again.
> delete(lon)
> delete(segments)
> delete(geometry)
> delete(segsDims)
> delete(geomDims)
> end do
>
>
> ;**********
> ; Vento *
> ;**********
>
> ; Risorse
>
>
> vento_direzione = plot_generale
>
>
>
> vento_direzione@vcRefLengthF = 0.06 ; define length of vec ref
> ;vento_direzione@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
> vento_direzione@vcRefAnnoArrowLineColor = "black" ; change ref
> vector color
> vento_direzione@vcRefAnnoArrowUseVecColor = False ; don't use vec
> color for ref
>
> ;vento_direzione@vcGlyphStyle = "CurlyVector" ; turn on curley vectors
> vento_direzione@vcLineArrowColor = "Black" ; change vector color
> vento_direzione@vcLineArrowThicknessF = 2.0 ; change vector thickness
>
> vento_direzione@vcLevelSelectionMode = "ManualLevels" ; set how
> manage contour line levels
> vento_direzione@vcMinLevelValF = 70 ; max level do draw contour lines
> vento_direzione@vcMaxLevelValF = 210 ; min level do draw contour lines
> vento_direzione@vcLevelSpacingF = 70
>
>
> vento_velocita = plot_generale
>
>
>
>
>
> NTIMES = dimsizes(gfs->time) ; number of times in the file
> U = gfs->ugrd_4572m * 3.6
> V = gfs->vgrd_4572m * 3.6 ; Convective rain
>
> VELOCITA_VENTO =
> round(((gfs->vgrd_4572m)/sin(atan2((gfs->vgrd_4572m),(gfs->ugrd_4572m))))*3.6,3)
>
>
> do it = 0,NTIMES-1 ;NTIMES-1
> testo_titolo_plot =
> gsn_create_text(wks,"Velocita~H-13V2F35~A~FV-2H3~ e direzione del
> vento",testo_titolo)
> draw(testo_titolo_plot)
> testo_modello_plot = gsn_create_text(wks,(run),testo_modello)
> draw(testo_modello_plot)
> testo_data_plot = gsn_create_text(wks,
> "Validita~H-13V2F35~A~FV-2H3~: "
> +(Nuova_data(ut_string(gfs->time(it),"%Y%N%D%H"))),testo_data)
> draw(testo_data_plot) ; Draw text box.
> precipitazioni_plot =
> gsn_csm_vector(wks,U(it,:,:),V(it,:,:),vento_direzione)
> overlay(map,precipitazioni_plot)
> draw(map)
> vento_velocita_plot =
> gsn_csm_contour(wks,VELOCITA_VENTO(it,:,:),vento_velocita)
> overlay(map,vento_velocita_plot)
> draw(map)
>
>
>
> destroy(precipitazioni_plot)
>
> frame(wks)
> end do
> do it = 0,NTIMES-1 ;NTIMES-1
>
>
> if (it.le.8) then
> system("mv Temperatura_2m_file_temporaneo.00000" + (it+1) + ".png" +
> " Temperatura0" + (it) + "_" + (run_nome_file) + ".png")
> end if
> if (it.eq.9) then
> system("mv Temperatura_2m_file_temporaneo.0000" + (it+1) + ".png" +
> " Temperatura0" + (it) + "_" + (run_nome_file) + ".png")
> end if
> if (it.ge.10) then
> system("mv Temperatura_2m_file_temporaneo.0000" + (it+1) + ".png" +
> " Temperatura" + (it) + "_" + (run_nome_file) + ".png")
> end if
> end do
> end
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk