Wind speed problem with no wind

From: ugo merlini <ugomerlini_at_nyahnyahspammersnyahnyah>
Date: Sun Oct 23 2011 - 11:00:20 MDT

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
Received on Sun Oct 23 11:00:30 2011

This archive was generated by hypermail 2.1.8 : Mon Oct 24 2011 - 09:29:35 MDT