Dear NCL,
I am trying to fix the issue with the joined Daymet grid, I used to
following script to generate a new file with the corrected coordinates:
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
Variable = "prcp"
nx = 1551
ny = 1872
N = nx * ny
;full_coordinate.dat
file1="full_coordinate.dat" ;the coordinate file
z3 = asciiread(file1,(/N,4/),"double")
lat1d=z3(:,2)
lon1d=z3(:,3)
dim2=(/ny,nx/)
lat2d=onedtond(lat1d,dim2)
lon2d=onedtond(lon1d,dim2)
printVarSummary(lat2d)
printVarSummary(lon2d)
;--find the zero value in the 2d lat and lon from the merged grid file
;NE_all_1980_prcp.nc
file2="NE_all_1980_" + Variable + ".nc"
ofile= "NE_all_1980_correct_" + Variable + ".nc"
o = addfile(ofile, "c")
f = addfile(file2,"r")
lat = f->lat
lon = f->lon
p = f->prcp
printVarSummary(p)
printVarSummary(lat)
printVarSummary(lon)
;--fill the zero with the generated coordinates
;If lat and lon are the same size as lat2d and lon2d, then you can use
"where", instead of a do loop:
lat_new = where(lat.eq.0,lat2d,lat)
lon_new = where(lon.eq.0,lon2d,lon)
o->prcp=f->prcp
;o->x=f->x
;o->y=f->y
o->lat=lat_new
o->lon=lon_new
end
The errors are:
ncendef: ncid 65536: NetCDF: One or more variable sizes violate format
constraints
fatal:NetCDF: Operation not allowed in define mode: error attempting to
write variable (time) to file
(/data/ecr/yangping/DAYMET/Grid/1KM/prcp/NE_all_1980_correct_prcp.nc)
warning:FileWriteVarVar: Could not write coordinate variable (time) to file
(NE_all_1980_correct_prcp), continuing anyway
ncredef: ncid 65536: NetCDF: Operation not allowed in define mode
ncendef: ncid 65536: NetCDF: One or more variable sizes violate format
constraints
fatal:NetCDF: Operation not allowed in define mode: error attempting to
write variable (y) to file
(/data/ecr/yangping/DAYMET/Grid/1KM/prcp/NE_all_1980_correct_prcp.nc)
warning:FileWriteVarVar: Could not write coordinate variable (y) to file
(NE_all_1980_correct_prcp), continuing anyway
ncredef: ncid 65536: NetCDF: Operation not allowed in define mode
ncendef: ncid 65536: NetCDF: One or more variable sizes violate format
constraints
fatal:NetCDF: Operation not allowed in define mode: error attempting to
write variable (x) to file
(/data/ecr/yangping/DAYMET/Grid/1KM/prcp/NE_all_1980_correct_prcp.nc)
warning:FileWriteVarVar: Could not write coordinate variable (x) to file
(NE_all_1980_correct_prcp), continuing anyway
warning:["Execute.c":7743]:Execute: Error occurred at or near line 48 in
file assign_coordinate.ncl
I don't understand why I can not do this? Is there a way that I can write
back the correct coordinate to the original file(without change other
value)?
Looking forward to hearing from you.
Best Regards,
Ping
On Mon, Nov 5, 2012 at 7:51 PM, Dave Allured <dave.allured@noaa.gov> wrote:
> 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
>
-- Ping Yang, Ph.D. CUNY Environmental Crossroads Initiative Marshak Science Building - Suite 925 The City College of New York - CUNY 160 Convent Avenue, New York NY 10031 Phone: 212-650-5769 Fax: 212-650-7064
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 19 14:00:36 2012
This archive was generated by hypermail 2.1.8 : Wed Nov 21 2012 - 11:16:05 MST