Re: About the new ESMF regridding !

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Mon Nov 05 2012 - 09:47:30 MST

Dear Echo,

The warning you got can be safely ignored. It's just telling you that the destination grid you have selected is slightly larger than your original source grid.

This would happen, for example, if your original grid went from -30 to 30 latitude, and 0 to 90 longitude, and then your new grid went from -40 to 40 latitude and 0 to 90 longitude.

I think your problem is with the use of tas_regrid and tas_new. You have:

 tas_regrid = ESMF_regrid_with_weights(tas,wgtFile,Opt)
. . .
    tas_new = tas_regrid

and then you do this:

 tas_new = new((/ntim,nlat,mlon/),float,UNDEF)

By doing the "new" command, you have effectively overwritten all the values that were in tas_new. Remove that line, and I think you should be okay.

Why do you have "tas_new" anyway? Can't you use tas_regrid directly?

Also, why are you making yet another copy of "tas_regrid" with:

  TAS = new((/ntim,nlat,mlon/),float,UNDEF)
  TAS(:,:,:) = tas_regrid(:,:,:)

You now have two copies of the same array, which can take up quite a bit more memory, and possibly time if your arrays are big.

 If you are using NCL V6.1.0 (not V6.1.0-beta), then "tas_regrid" should automatically have all the coordinate arrays attached, so then you won't need all of the code for assigning metadata.

I think your script should look something like this (UNTESTED):

 
begin
;---Input file
diri = "./"
fili = "ACCESS1-0/Historial/r1i1p1/tas_Amon_ACCESS1-0_historical_r1i1p1_185001-200512.nc"

 
    year = 1901 ; TDEF
    ntim = 1260 ; time steps
    nmos = 12
    UNDEF = 1E+20
 
    srcFileName = diri+fili
 
;---Output (and input) files
    srcGridName = "ACCESS1-0_SCRIP.nc" ;
    dstGridName = "World1deg_SCRIP.nc"
    wgtFile = "ACCESS1-0_2_World_SCRIP.nc"
 
;---Retrie ve data from BCC Grid
    sfile = addfile(srcFileName,"r")
 
    nlat = 180 ; YDEF
    mlon = 360 ; XDEF
 
  tas = sfile->tas(612:1871,:,:) ;
;---Set to True if you want to skip any of these steps
 
    SKIP_BCC_SCRIP_GEN = False
    SKIP_WORLD_SCRIP_GEN = False
    SKIP_WGT_GEN = False
 
;----------------------------------------------------------------------
; Convert BCC to SCRIP file.
;----------------------------------------------------------------------
    nav_lat = sfile->lat
    nav_lon = sfile->lon
 if(.not.SKIP_BCC_SCRIP_GEN) then
;---Convert to an SCRIP Convention file.
      Opt = True
      Opt@ForceOv erwrite = True
      Opt@PrintTimings = True
      Opt@Mask2D = where(.not.ismissing(tas(0,:,:)),1,0)
 
      rectilinear_to_SCRIP(srcGridName,nav_lat,nav_lon,Opt)
 
;---Clean up
      delete(Opt)
    end if
 
;----------------------------------------------------------------------
; Convert 1 degree world grid to SCRIP file
;----------------------------------------------------------------------
    if(.not.SKIP_WORLD_SCRIP_GEN)
      Opt = True
      Opt@LLCorner = (/-89.75d, 0.00d /)
      Opt@URCorner = (/ 89.75d, 359.75d /)
      Opt@ForceOverwrite = True
      Opt@PrintTimings = True
      Opt@Title = "World grid 1x1 degree resolution"
 
      latlon_to_SCRIP(dstGridName,"1deg",Opt)
 
;---Clean up
      delete(Opt)
    end if
 
;----------------------------------------------------------------------
; Generate interpolation weights for BCC Grid to World Grid
;----------------------------------------------------------------------
  &nbs p; if(.not.SKIP_WGT_GEN) then
      Opt = True
      Opt@SrcESMF = False
      Opt@DstESMF = False
  Opt@ForceOverwrite = True
      Opt@PrintTimings = True
 
      ESMF_regrid_gen_weights(srcGridName, dstGridName, wgtFile, Opt)
 
;---Clean up
      delete(Opt)
    end if
 
;----------------------------------------------------------------------
; Interpolate data from BCC to World 1-degree grid.
;----------------------------------------------------------------------
 
    Opt = True
; Opt@Debug = True
    Opt@PrintTimings = True
 
    tas_regrid = ESMF_regrid_with_weights(tas,wgtFile,Opt)
    printVarSummary(tas_regrid)
 
;---Clean up
    delete(Opt)
 
 
;------output regridding data-------
  NC = True
 
    year = 1901 ; ; TDEF
    ntim = 1260 ; time steps
    nmos = 12
    UNDEF = 1E+20
    nlat = 180 ; YDEF
    mlon = 360 ; XDEF

    yyyy = ispan(1901,2005,1)
    mm = ispan(1,12,1)
    ny = dimsizes(yyyy)
    nm = dimsizes(mm)
    years = conform_dims((/ny,nm/),yyyy,0)
    months = conform_dims((/ny,nm/),mm,1)
    time = sprinti("%4i",ndtooned(years)) + \
             sprinti("%0.2i", ndtooned(months))
 
  time!0 = "time"
  time@long_name = "time"
  time@units = "yyyymm"
  time&time = time
 
; delete(tas_regrid&time) ; might need this
  tas_regrid&time = time ; Assign new time values.
  tas_regrid@long_name = "Near-Surface Temperature" ; Vars
  tas_regrid@units = "K"
  tas_regrid@_FillValue = 1E+20 ; UNDEF
 
  printVarSummary(tas_regrid)
  print ("min(tas_regrid)="+min(tas_regrid))
  print ("max(tas_regrid)="+max(tas_regrid))
 
  if (NC) then
     diro = "./" ; Output directory
     filo = "tas_Amon_ACCESS1-0_historical_r1i1p1_185001-201212-1x1.nc"
     system("/bin/rm -f " + diro + filo) ; remove if exists
 
     ncdf = addfile (diro+filo, "c")
 
; make time and UNLIMITED dimension ; recommended for most applications
 
        filedimdef(ncdf,"time",-1,True)
 
; output vari ables directly
 
       ncdf->TAS = tas_regrid
 
  end if
exit
 
end

--Mary

On Nov 3, 2012, at 3:45 AM, EchoChen wrote:

> Dear NCL community,
>
>
> First, I want to thank MARY for her suggestion. I have tried Mary's ways ,but I failed in outputing data. And I don't know how to deal with some errors(WARNING).
> Scripts and results below are in MARY's way and in an other way with a loop . Kindly hope someone give me more suggestions.
> Thanks for your attention.
>
> Echo
>
>
> NO.1(Mary's way):
>
> begin
> ;---Input file
> diri = "./"
> fili = "ACCESS1-0/Historial/r1i1p1/tas_Amon_ACCESS1-0_historical_r1i1p1_185001-200512.nc"
>
>
> year = 1901 ; TDEF
> ntim = 1260 ; time steps
> nmos = 12
> UNDEF = 1E+20
>
> srcFileName = diri+fili
>
> ;---Output (and input) files
> srcGridName = "ACCESS1-0_SCRIP.nc" ;
> dstGridName = "World1deg_SCRIP.nc"
> wgtFile = "ACCESS1-0_2_World_SCRIP.nc"
>
> ;---Retrie ve data from BCC Grid
> sfile = addfile(srcFileName,"r")
>
> nlat = 180 ; YDEF
> mlon = 360 ; XDEF
>
> ; tas_new = new((/ntim,nlat,mlon/),float,UNDEF)
> ; tos =new((/ntim,nlat,mlon/),float,UNDEF)
>
> tas = sfile->tas(612:1871,:,:) ;
> ;---Set to True if you want to skip any of these steps
>
> SKIP_BCC_SCRIP_GEN = False
> SKIP_WORLD_SCRIP_GEN = False
> SKIP_WGT_GEN = False
>
> ;----------------------------------------------------------------------
> ; Convert BCC to SCRIP file.
> ;----------------------------------------------------------------------
> nav_lat = sfile->lat
> nav_lon = sfile->lon
> if(.not.SKIP_BCC_SCRIP_GEN) then
> ;---Convert to an SCRIP Convention file.
> Opt = True
> Opt@ForceOv erwrite = True
> Opt@PrintTimings = True
> Opt@Mask2D = where(.not.ismissing(tas(0,:,:)),1,0)
>
> rectilinear_to_SCRIP(srcGridName,nav_lat,nav_lon,Opt)
>
> ;---Clean up
> delete(Opt)
> end if
>
> ;----------------------------------------------------------------------
> ; Convert 1 degree world grid to SCRIP file
> ;----------------------------------------------------------------------
> if(.not.SKIP_WORLD_SCRIP_GEN)
> Opt = True
> Opt@LLCorner = (/-89.75d, 0.00d /)
> Opt@URCorner = (/ 89.75d, 359.75d /)
> Opt@ForceOverwrite = True
> Opt@PrintTimings = True
> Opt@Title = "World grid 1x1 degree resolution"
>
> latlon_to_SCRIP(dstGridName,"1deg",Opt)
>
> ;---Clean up
> delete(Opt)
> end if
>
> ;----------------------------------------------------------------------
> ; Generate interpolation weights for BCC Grid to World Grid
> ;----------------------------------------------------------------------
> &nbs p; if(.not.SKIP_WGT_GEN) then
> Opt = True
> Opt@SrcESMF = False
> Opt@DstESMF = False
> Opt@ForceOverwrite = True
> Opt@PrintTimings = True
>
> ESMF_regrid_gen_weights(srcGridName, dstGridName, wgtFile, Opt)
>
> ;---Clean up
> delete(Opt)
> end if
>
> ;----------------------------------------------------------------------
> ; Interpolate data from BCC to World 1-degree grid.
> ;----------------------------------------------------------------------
>
> Opt = True
> ; Opt@Debug = True
> Opt@PrintTimings = True
>
> tas_regrid = ESMF_regrid_with_weights(tas,wgtFile,Opt)
> printVarSummary(tas_regrid)
>
> tas_new = tas_regrid
> printVarSummary(tas_new)
>
> ;---Clean up
> delete(Opt)
>
>
> ;------output regridding data-------
> NC = True
>
> year = 1901 ; ; TDEF
> ntim = 1260 ; time steps
> nmos = 12
> UNDEF = 1E+20
> nlat = 180 ; YDEF
> mlon = 360 ; XDEF
>
> tas_new = new((/ntim,nlat,mlon/),float,UNDEF)
> ;==================== =========
> ; create coordinate variables
> ;============================
> lon=fspan(0.75,359.75,360)
> lon!0 = "lon"
> lon@long_name = "longitude"
> lon@units = "degrees_east"
> lon&lon = lon
>
> lat=fspan(-89.75,89.75,180)
> lat!0 = "lat"
> lat@long_name = "latitude"
> lat@units = "degrees_north"
> lat&lat = lat
>
> years = sprinti("%4i",ispan(1901,2005,1)) ; write integers
> months = sprinti("%0.2i", ispan(1,12,1)) ; ditto for
>
> time = new(dimsizes(years)*dimsizes(months),"integer") ; create
> k=0
> do i=0,dimsizes(years)-1
> do j = 0, dimsizes(months)-1
> time(k)=stringtointeger(years(i) +months(j)) ; combine to
> k=k+1
> end do
> end do
>
> time!0 = "time"
> time@long_name = "time"
> time@units = "yyyymm"
> time&time = time
>
> tas_new!0 = "time"
> tas_new!1 = "lat"
> tas_new!2 = "lon"
>
> tas_new&time = time
> tas_new&lat = lat
> tas_new&lon = lon
>
> tas_new@long_name = "Near-Surface Temperature" ; Vars
> tas_new@units = "K"
> tas_new@_FillValue = 1E+20 ; UNDEF
>
> print VarSummary(tas_new)
> print ("min(tas_new)="+min(tas_new))
> print ("max(tosregrid)="+max(tas_new))
>
> TAS = new((/ntim,nlat,mlon/),float,UNDEF)
> TAS(:,:,:) = tas_new(:,:,:)
> copy_VarMeta(tas_new,TAS)
>
> if (NC) then
> diro = "./" ; Output directory
> filo = "tas_Amon_ACCESS1-0_historical_r1i1p1_185001-201212-1x1.nc"
> system("/bin/rm -f " + diro + filo) ; remove if exists
>
> ncdf = addfile (diro+filo, "c")
>
> ; make time and UNLIMITED dimension ; recommended for most applications
>
> filedimdef(ncdf,"time",-1,True)
>
> ; output vari ables directly
>
> ncdf->TAS = TAS
>
> end if
> exit
>
> end
>
>
> (0) =====> CPU Elapsed Time: rectilinear_to_SCRIP: 3.63145 seconds <=====
> (0) =====> CPU Elapsed Time: latlon_to_SCRIP: 8.18476 seconds <=====
> (0) =====> CPU Elapsed Time: ESMF_regrid_gen_weights: 0.299954 seconds <=====
> (0) ESMF_regrid_with_weights: warning: destination grid is not completely
> (0) covered by the source grid.
> (0) =====> CPU Elapsed Time: ESMF_regrid_with_weights: 14.9517 seconds <=====
>
>
> The result is that there is no value in the new ......1x1.nc file except lat/lon information. I don't know why that warning: destination grid is not completely covered by the source grid. Actually, I don't know how to change Opt@Mask2D and Opt@LLCorner Opt@URCorner .
> I change Opt@LLCorner to (/-90d, 0.00d /) and Opt@URCorner to (/ 90d, 360d /) but it can't work. Why? Can I use the Opt@Mask2D to mask ocean data or else?
>
>
>
> NO.2(loop)
>
> begin
> ;---Input file
> diri = "./"
> fili = "tas_Amon_ACCESS1-0_historical_r1i1p1_185001-200512.nc"
>
>
> year = 1901 ; TDEF
> ntim = 1260 ; time steps
> nmos = 12
> UNDEF = 1E+20
> srcFileName = diri+fili
> ;---Output (and input) files
> srcGridName = "ACCESS1-0_SCRIP.nc" ;
> dstGridName = "World1deg_SCRIP.nc"
> wgtFile = "ACCESS1-0_2_World_SCRIP.nc"
> ;---Retrieve data from BCC Grid
> sfile = addfile(srcFileName,"r")
> nlat = 180 ; YDEF
> mlon = 360 ; XDEF
> tas_new = new((/ntim,nlat,mlon/),float,UNDEF)
>
> do nt = 0,ntim-1
> tos = sfile->tas(nt+612,:,:) ;624 change to 612
> ;---Set to True if you want to skip any of these steps
> SKIP_BCC_SCRIP_GEN = False
> SKIP_WORLD_SCRIP_GEN = False
> SKIP_WGT_GEN = False
>
> ;----------------------------------------------------------------------
> ; Convert BCC to SCRIP file.
> ;----------------------------------------------------------------------
> nav_lat = sfile->lat
> nav_lon = sfile->lon
> if(.not.SKIP_BCC_SCRIP_GEN) then
> ;---Convert to an SCRIP Convention file.
> Opt = True
> Opt@ForceOverwrite = True
> Opt@PrintTimings = True
> Opt@Mask2D = where(.not.ismissing(tos),1,0)
> rectilinear_to_SCRIP(srcGridName,nav_lat,nav_lon,Opt)
> ;---Clean up
> delete(Opt)
> end if
> ;----------------------------------------------------------------------
> ; Convert 1 degree world grid to SCRIP file
> ;----------------------------------------------------------------------
> if(.not.SKIP_WORLD_SCRIP_GEN)
> Opt = True
> Opt@LLCorner = (/-89.75d, 0.00d /)
> Opt@URCorner = (/ 89.75d, 359.75d /)
> Opt@ForceOverwrite = True
> Opt@PrintTimings = True
> Opt@Title &nbs p; = "World grid 1x1 degree resolution"
> latlon_to_SCRIP(dstGridName,"1deg",Opt)
> ;---Clean up
> delete(Opt)
> end if
> ;----------------------------------------------------------------------
> ; Generate interpolation weights for BCC Grid to World Grid
> ;----------------------------------------------------------------------
> if(.not.SKIP_WGT_GEN) then
> Opt = True
> Opt@SrcESMF = False
> Opt@DstESMF = False
> Opt@ForceOverwrite = True
> Opt@Pri ntTimings = True
> ESMF_regrid_gen_weights(srcGridName, dstGridName, wgtFile, Opt)
> ;---Clean up
> delete(Opt)
> end if
> ;----------------------------------------------------------------------
> ; Interpolate data from BCC to World 1-degree grid.
> ;----------------------------------------------------------------------
> Opt = True
> ; Opt@Debug = True
> Opt@PrintTimings = True
> tos_regrid = ESMF_regrid_with_weights(tos,wgtFile,Opt)
> printVarSummary(tos_regrid)
> tas_new(nt,:,:) = tos_regrid(:,:)
> printVarSummary(tas_new)
> ;---Clean up
> delete(Opt)
>
> end do
> ;------output regridding data-------
> NC = True
> ;=============================
> ; create coordinate variables
> ;============================
> lon=fspan(0.75,359.75,360)
> lon!0 = "lon"
> lon@long_name = "longitude"
> lon@units = "degrees_east"
> lon&lon = lon
> lat=fspan(-89.75,89.75,180)
> lat!0 = "lat"
> lat@long_name = "latitude"
> lat@units = "degrees_north"
> lat&lat = lat
> years = sprinti("%4i",ispan(1901,2005,1)) ; write integers
> months = sprinti("%0.2i", ispan(1,12,1)) ; ditto for
> time = new(dimsizes(years)*dimsizes(months),"integer") ; create
> k=0
> do i=0,dimsizes(years)-1
> do j = 0, dimsizes(months)-1
> time(k)=stringtointeger(years(i)+months(j)) ; combine to
> k=k+1
> end do
> end do
> time!0 = "time"
> time@long_name = "time"
> time@units = "yyyymm"
> time&time = time
> tas_new!0 = "time"
> tas_new!1 = "lat"
> tas_new!2 = "lon"
> tas_new&time = time
> tas_new&lat = lat
> tas_new&lon = lon
> tas_new@long_name = "Near-Surface Temperature" ; Vars
> tas_new@units = "K"
> tas_new@_FillValue = 1E+20 ; UNDEF
> printVarSummary(tas_new)
> print ("min(tas_new)="+min(tas_new))
> print ("max(tosregrid)="+max(tas_new))
> TAS = new((/ntim,nlat,mlon/),float,UNDEF)
> TAS(:,:,:) = tas_new(:,:,:)
> copy_VarMeta(tas_new,TAS)
> if (NC) then
> diro = "./" ; Output directory
> filo = "tas_Amon_ACCESS1-0_historical_r1i1p1_185001-201212-1x1.nc"
> system("/bin/rm -f " + diro + filo) ; remove if exists
> ncdf = addfile (diro+filo, "c")
> ; make time and UNLIMITED dimension ; recommended for most applications
> filedimdef(ncdf,"time",-1,True)
> ; output variables directly
> ncdf->TAS = TAS
> end if
> exit
>
>
> Result: It seems to be correct but takes a long time (more than 3 hours).
>
>
>
>
>
>
>
>
> > Subject: Re: [ncl-talk] About the new ESMF regridding !
> > From: haley@ucar.edu
> > Date: Fri, 2 Nov 2012 09:40:51 -0600
> > CC: ncl-talk@ucar.edu
> > To: chenxiaochen_ams@hotmail.com
> >
> > Dear Echo,
> >
> > Welcome to the NCL community.
> >
> > In the ESMF_all_7.ncl script, the "test_CNRM.nc" file actually only has one time step, so that's partially why it's only reading the one time step.
> >
> > However, if this file had multiple time steps, then you can easily change to multiple time steps by changing this line:
> >
> > tos = sfile->tos(0,:,:)
> >
> > to:
> >
> > tos = sfile->tos
> >
> >
> > Assuming the mask is the same for each timestep, then you also need to change this line:
> >
> > Opt@Mask2D = where(.not.ismissing(tos),1,0)
> >
> > to:
> >
> > Opt@Mask2D = where(.not.ismissing(tos(0,:,:)),1,0)
> > ;
> > When you get to this call:
> >
> > tos_regrid = ESMF_regrid_with_weights(tos,wgtFile,Opt)
> >
> > It should automatically regrid across however many dimensions you have. "tos_regrid" will have the same leftmost dimensions as "tos", and the rightmost dimensions will be whatever grid you regridded it to.
> >
> > The important thing to note is that this method will only work if your grids are the same across all leftmost dimensions. If you have a source or destination grid that changes with each time or level or any other leftmost dimension, then you need to generate a separate weights file for each of these.
> >
> > If you look at example ESMF_all_6.ncl, which regrids data from a CMIP5 grid to a 1x1 grid, the same rules apply. This file has multiple timesteps and levels, so you would simply change two lines:
> >
> > thetao = sfile->thetao(0,0,:,:)
> >
> > to:
> >
> > thetao = sfile->thetao
> >
> > and:
> >
> > Opt@Mask2D = where(.not.ismissing(thetao),1,0)
> >
> > to:
> >
> > Opt@Mask2D = where(.not.ismissing(thetao(0,0,:,:)),1,0)
> >
> > In this case, then, since "thetao" is dimensioned 600 x 19 x 220 x 256 (time x lev x lat x lon), then "thetao_regrid" will be dimensioned 600 x 19 x 180 x 360, since 180 x 360 is the size of the destination (1x1) grid.
> >
> > Hope this helps clear things up.
> >
> > --Mary
> >
> >
> > On Nov 1, 2012, at 8:03 PM, ³ÂÏþ³¿ wrote:
> >
> > > Dear ncl communit,
> > >
> > >
> > > I am a new leaner. It's first time for me to send a email to the community. English is not my mother tougue, so maybe there are some mistakes in my email. Thank you for your understanding.
> > >
> > > I run the new ESMF regridding script- ESMF_all_7.ncl , and I succeed in plotting after I modify some instrctions.( http://www.ncl.ucar.edu/Ap plications/ESMF.shtml) I also scan other scripts in the ESMF regridding example pages for reference. But I want to regrid a 3D or 4D (or higher) data like A(:,:,:) and output it ,but all scripts just regrid a 2D data like A(0,:,:) or A(0,0,:,:).
> > >
> > > I have a large number of datum to regrid. My datum on rectilinear are from CMIP. The type of interpolation methods is bilinear.Do you have some suggestions for me? Maybe I should use cycling but I do not know how.
> > >
> > > Thanks for your attention!
> > >
> > >
> > >
> > > Echo
> > > _______________________________________________
> > > 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
Received on Mon Nov 5 09:47:40 2012

This archive was generated by hypermail 2.1.8 : Tue Nov 06 2012 - 15:05:49 MST