Re: Regridding Daymet Joined Tiles with ESMF_Regridding on NCL 6.1.0

Hello Ping, Mary, et al.

I just put in a request for fully populated reference grids of 2-D
coordinates with Daymet/ORNL. If they can provide, I will be able to
upgrade the original Daymet join script to make merged files without
gaps in the 2-D coordinates. This should fix these problems with ESMF
regridding and other NCL regridders that you might want to use.

I recommend this approach rather than regridding individual tiles.
The concept of the join script carries over to 1-D coordinates, but
there will be significant re-coding for several things that change.

If you wish to go ahead with 1-D coordinates, then mainly you would
remove everything concerned with variables "olat" and "olon" which are
the merged 2-D coordinates, no longer used. Determine the start and
end indices for overlays from the 1-D coordinates f1->lon, f1->lat,
etc., the same way that the original version used f1->x, f1->y, etc.
The use of overlay subscripts i1a, i1b, etc. should remain about the
same. Also review all of the output metadata, remove things not
needed, and add fixes as needed. The original script was very
specific for the layout and metadata of Daymet files.

--Dave

On Mon, Nov 5, 2012 at 3:58 PM, Ping Yang <pyang@ccny.cuny.edu> wrote:
> 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
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 5 17:51:58 2012