Re: Regridding Daymet Joined Tiles with ESMF_Regridding on NCL 6.1.0

Hi NCL,

I think there is something wrong with the source grid, as the error is that
> the grid corners were not correct. You can see this by looking at the
> "PET0.RegridWeightGen.Log" file:
>
> cat PET0.RegridWeightGen.Log
> 20121105 115739.811 INFO PET0 Running with ESMF Version 5.2.0rp2
> 20121105 115740.616 ERROR PET0 ESMF_Grid.F90:5238
> ESMF_GridCreateFrmScripDistGrd Wrong argument specified - - Bad corner
> array in SCRIP file
> 20121105 115740.616 ERROR PET0 ESMF_Grid.F90:5585
> ESMF_GridCreateFrmScripReg Wrong argument specified - Internal subroutine
> call returned Error
> 20121105 115810.964 INFO PET0 Running with ESMF Version 5.2.0rp2
>
> When I try to change the interpolation method to "bilinear", it just runs
> forever.
>
> BTW, it's a good idea to set:
>
> Opt@Debug = True
>
Good suggestion.

>
> when you are initially writing the regridding script, so you can get some
> more debug information printed out.
>
> How was the source grid created? Is this one tile, or several tiles put
> together in some fashion?
>
> If it's the latter, then this is likely the source of your problem.
>
The source grid is created by joining several tiles together, I think
that possibly the problem come from.

I am thinking a way to avoid this problem, eventually I need are lat/lon
rectilinear grids for our model input. The way I am thinking is: regridding
the curvilinear grid (daymet tiles) into rectilinear grid first(or just
reprojecting it) and then use a script to join them.

Is there an existing example for joining two adjacent rectilinear grid
files? I am trying to follow the daymet joining method but I got lost.

Following is my script for joining two adjacent regridded daymet files:

   f1 = addfile (file1, "r") ; open input files
   f2 = addfile (file2, "r")

   x1 = f1->lon ; read 1-D coordinate vectors
   y1 = f1->lat
   x2 = f2->lon
   y2 = f2->lat

   print(x1)
   print(x2)
   print(y1)
   print(y2)
   ;ind--return the indices where the input is True
   i1a = ind (x2 .eq. x1(0)) ; compute relative offsets
   i2a = ind (x1 .eq. x2(0)) ; between the two input grids
   j1a = ind (y2 .eq. y1(0))
   j2a = ind (y1 .eq. y2(0))

; print(i1a)
; print(i2a)
; print(j1a)
; print(j2a)

   if ( all (ismissing ((/ i1a, i2a /))) \
    .or. all (ismissing ((/ j1a, j2a /))) ) then
      print ("*** No common X/Y grid point, join fails. Abort.")
      print ("*** File 1 = " + file1)
      print ("*** File 2 = " + file2)
      status_exit (1)
   end if

   i1a = max ((/ 0, i1a /)) ; compute starting offsets
   i2a = max ((/ 0, i2a /)) ; with missing value trick
   j1a = max ((/ 0, j1a /))
   j2a = max ((/ 0, j2a /))

   print(i1a)
   print(i2a)
   print(j1a)
   print(j2a)

   print ("i1a, i2a =" + sprinti ("%6i", i1a) + sprinti ("%6i", i2a))
   print ("j1a, j2a =" + sprinti ("%6i", j1a) + sprinti ("%6i", j2a))

   i1b = i1a + dimsizes (x1) - 1 ; compute ending offsets
   i2b = i2a + dimsizes (x2) - 1
   j1b = j1a + dimsizes (y1) - 1
   j2b = j2a + dimsizes (y2) - 1

   print ("i1b, i2b =" + sprinti ("%6i", i1b) + sprinti ("%6i", i2b))
   print ("j1b, j2b =" + sprinti ("%6i", j1b) + sprinti ("%6i", j2b))

   nxout = 1 + max ((/ i1b, i2b /)) ; combined grid size
   nyout = 1 + max ((/ j1b, j2b /))

   print ("nxout, nyout = " + nxout + ", " + nyout)

   ntimes = dimsizes (f1->time) ; allocate output arrays
; dims2 = (/ nyout, nxout /)
   dims3 = (/ ntimes, nyout, nxout /)

   sample = f1->$var$(0,0,0)

    olat = new (nyout, typeof (f1->lat),"No_FillValue")
    olon = new (nxout, typeof (f1->lon),"No_FillValue")

; olat = new (dims2, typeof (f1->lat), "No_FillValue")
; olon = new (dims2, typeof (f1->lon), "No_FillValue")
   xout = new (dims3, typeof (sample), sample@_FillValue)
   printVarSummary(olat)
   printVarSummary(olon)

   ;olat(:,:) = 0 ; back fill 2-D coordinates with
   ;olon(:,:) = 0 ; zero as pseudo-missing value

   print ("Overlay first input grid into arrays.")

; olat(j1a:j1b,i1a:i1b) = f1->lat ; overlay 2-D coordinates
; olon(j1a:j1b,i1a:i1b) = f1->lon
   olat(i1a:i1b) = f1->lat
   olon(i1a:i1b) = f2->lon

   xout(:,i1a:i1b,i1a:i1b) = f1->$var$ ; overlay main array

   print ("Overlay second input grid into arrays.")

  ;olat(j2a:j2b,i2a:i2b) = (/ where ((f2->lat .eq. 0), \
    ; olat(j2a:j2b,i2a:i2b), f2->lat) /) ; overlay 2-D coordinates; use
   ; mask to preserve existing coords
   olat(i2a:i2b) = (/where((f2->lat .eq. 0, olat(i2a:i2b),f2->lat)/)

   ;olon(j2a:j2b,i2a:i2b) = (/ where ((f2->lon .eq. 0), \
     ; olon(j2a:j2b,i2a:i2b), f2->lon) /)
   olon(i2a:i2b) = (/where((f2->lon .eq. 0, olon(i2a:i2b), f2->lon)/)
   ;xout(:,j2a:j2b,i2a:i2b) = (/ where (ismissing (f2->$var$), \
    ; xout(:,j2a:j2b,i2a:i2b), f2->$var$) /)
      ; overlay main array; use mask
; to preserve existing data

   xout(:,i2a:i2b,j2a:j2b) =
(/where(ismissing(f2->$var$),xout(:,i2a:i2b,j2a:j2b),f2->$var$)/)

   print ("Write output file: " + outfile)

   if (isfilepresent (outfile)) then ; overwrite any previous file
      system ("rm " + outfile)
   end if

   out = addfile (outfile, "c") ; create new output file

I will be appreciated if a specialist can point out the problem. Actually I
didn't figure out how to deal with the overlaying area, e.g. how to assign
the lon and lat value for it.

Many thanks,

Ping
Received on Mon Nov 5 15:59:41 2012