Re: Regridding question

From: Dennis Shea (shea AT cgd.ucar.edu)
Date: Fri Oct 07 2005 - 15:59:56 MDT

  • Next message: S. Yi: "problem of using gsn_attach_plots"

    Hello,

    Pls see attached text. It will do the interpolation.

    A few comments for a person "new to NCL"

    [1] pc(lat_2,lon_2) ; [ 90N -> 90S]

        In NCL, use array syntax to reverse the order of
        a dimension.
        
        pc = pc(::-1,:) ; [ 90S -> 90N] done "in place"
        
    [2] There are *negative* prc values.
        This may indicate that values originated with
        a spectral model.
        
    [3] Become familiar with
           [a] contributed.ncl
           [b] the "Application Examples"
           
        They are your friends :-)
        
        
    [4] GRIB files generqally do not come with the ".grb"
        extension. This is not required for the file name.
        If the file name is
        
             gfs.t00z.pgrbf192
             
        Then you can use:
        
         f = addfile("gfs.t00z.pgrbf192.grb","r")
         
         NCL will look for the file " gfs.t00z.pgrbf192.grb",
         If this is not found, NCL will look for the file
         "gfs.t00z.pgrbf192" but 'remember' to treat it
         as a GrIB file.
         
         This goes for netCDF, hdf, hdfeos, etc
         Onny in addfile must the extension occure. This
         tells NCL how to read the file
         
        
    good luck
    D
         
    ===========================================================
    Hi Guys,
    I am new to NCL and am trying to convert a 73x144 lat lon grid to a 181x360 lat
    lon grid using the linint2 function. Heres what I have done so far:

    gfile=addfile("gfs.t00z.pgrbf192.grb","r") ;read a GFS model grib file
    load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
    load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    pc=gfile->A_PCP_2_SFC_acc12h ; read the 12 hour Total surface ppt (lat| 73, lon|
    144)

    ; Since linint2 expects lat array in ascending order hence:
    pcflip = new((/73,144/),float)
    k=72
    do i = 0,72
    pcflip(i,:) = pc(k,:)
       k = k-1
    end do
    newlon=fspan(0.,359,360) ; define new longitude
    newlat=fspan(-90.,90,181) ; define new output latitude
    pc_regrid = linint2_Wrap(pcflip&lon_2,pcflip&lat_2,pcflip,False,newlon,newlat,0)

     ; Then flip the data back again along the latitude
     k = 180
    do i = 0,180
    pcregrid_flip(i,:) = pc_regrid(k,:)
      k = k-1
    end do

    So my final regridded varibale is pcregrid_flip
    I am not sure i am doing the right thing here? the plots of the original and
    regridded variables looks somewhat different which maybe due to interpolation
    but i just want to make sure from you guys. Also shuold i use f2sfh insetad of
    linint2??
    Any suggestion will be most welcome
    Thanks

    Sudipta Sarkar


    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"

    begin
      diri = "/project/cas/shea/GRIB/"
      fili = "gfs.t00z.pgrbf192"
      var = "A_PCP_2_SFC_acc12h"

      f = addfile(diri+fili+".grb","r")
      pc = f->$var$ ; North to South
     ;print(f) ; like ncdump [overview of file]

      printVarSummary(pc)
      printMinMax(pc, True) ; This has *negative* values
    ;************************************************
    ; bilinear interpolation
    ;************************************************
      nlat = 181
      mlon = 360
      latNew= latGlobeF (nlat, "lat", "latitude", "degrees_north")
      lonNew= lonGlobeF (mlon, "lon", "longitude", "degrees_east")

      pc = pc(::-1,:) ; South to North for linint2
      printVarSummary(pc) ; Note coordinate arrays
      printMinMax(pc, True)

      pcNew = linint2_Wrap(pc&lon_2,pc&lat_2,pc,True ,lonNew,latNew,0)

      printVarSummary(pcNew)
      printMinMax(pcNew, True)

      pcNew = pcNew(::-1,:) ; North to South
      printVarSummary(pcNew)

      itime = pc@initial_time
      ftime = pc@forecast_time
      
    ;************************************************
    ; create plots
    ;************************************************
      wks = gsn_open_wks("x11","ssarkar_prc")
     ;gsn_define_colormap(wks,"gui_default")

      colors = (/"white","black" \ ; back/fore ground
                ,"blue","azure1","beige","lavender" \
                ,"PaleGreen","SeaGreen3","yellow" \
                ,"magenta","HotPink","Red"/) ; choose colors for color map
      gsn_define_colormap(wks, colors) ; generate new color map

      plot = new(2,graphic) ; create a plot array
      
      res = True
      res@gsnDraw = False ; don't draw
      res@gsnFrame = False ; don't advance frame
      res@cnInfoLabelOn = False ; turn off cn info label
      res@cnFillOn = True ; turn on color
      res@lbLabelBarOn = False ; turn off individual cb's

      res@cnLevelSelectionMode= "ExplicitLevels" ; set explicit contour levels
                                                   ; specify unequal contour levels
      res@cnLevels = (/ 0.0, 1.0, 3.0, 5.0 ,10 ,15 ,25 ,50, 75 /)

      res@gsnCenterString = "2.5"
      plot(0) = gsn_csm_contour_map_ce(wks,pc,res)
      res@gsnCenterString = "1x1"
      plot(1) = gsn_csm_contour_map_ce(wks,pcNew,res)
    ;************************************************
    ; create panel
    ;************************************************
      resP = True ; modify the panel plot
      resP@txString = "itime="+itime+" ftime="+ftime
      resP@gsnPanelLabelBar = True ; add common colorbar
      resP@gsnMaximize = True ; if "ps" will make large
      gsn_panel(wks,plot,(/2,1/),resP) ; now draw as one plot
    end


    _______________________________________________
    ncl-talk mailing list
    ncl-talk@ucar.edu
    http://mailman.ucar.edu/mailman/listinfo/ncl-talk



    This archive was generated by hypermail 2b29 : Fri Oct 07 2005 - 17:58:12 MDT