[Thread Prev][Thread Next][Index]
Re: Comparing Two Curvilinear Variables
Hi,
It seems that the transformation of curvilinear data to regular
grid is of general interst. Earlier Prasad posted a mail about the same
(http://ferret.pmel.noaa.gov/Ferret/Mail_Archives/fu_2003/msg01271.html)
and we made some offline discussions. Infact the answer to this problem
was there in the Ferret mail archive
(http://www.ferret.noaa.gov/Ferret/Mail_Archives/fu_2003/msg00286.html)
as a workaround suggested for SAMPLEXY_CURV function by James Orr.
James Orr kindly provided me all the details necessary for the
transformation of a variable on curvilinear co-ordinates to regular/
rectilinear co-ordinates. I just made a simple example based on
objective_analysis_demo.jnl to describe the method. Fine details are
given in the jnl(curv2regular_demo.jnl) file atached.
Suppose the curvilinear data is of the form (from Prasad's mail)
dimensions:
nx = 316 ;
ny = 332 ;
variables:
float lon(ny, nx) ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
float lat(ny, nx) ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
float ice_con(ny, nx) ;
ice_con:units = "%"
ie. ny*nx lat points, ny*nx lon points, ny*nx data(ice_con) points.
The steps involved to get this in regular grid are :
1. get 2D data of lat(ny,nx) & lon(ny,nx) to 1D variables
of nx*ny points using xsequence
2. define regular/rectilinear axes on which the data
is needed
3. Now we have 1D lat (nx*ny points), 1D lon (nx*ny points)
and 2D(in the above case; only x & y dimesions)
variable (nx*ny points) ...or we can say
data in a "scattered" form. Regrid this data to
a "regular grid" using SCAT2GRIDGAUSS_XY or SCAT2GRIDLAPLACE_XY
functions available in ferret. We need to provide
other necessary arguments like XSCALE/YSCALE/OFFSET
to the function according to the accuracy/need.
4. perform nessessary calculations/plotting using the new
variable on regular grid and save the output to
a netcdf file for future use (if needed).
Daryl almost reached the answer ; the problem was with the 2D array of
latitude and longitude passed to the SCAT2GRIDGAUSS_XY function.
I would like to thank James Orr for his quick response and detailed
reply. Hope James Orr will add his comments especially if I missed
something.
With Regards
Jaison
On Tue, 29 Jun 2004, Daryl Herzmann wrote:
> Hi Ferret!
>
> I am in a dill of a pickle and have been trying to think like a Ferret
> without much luck. Perhaps a kind person could help me out with a
> documentation reference I have probably missed :(
>
> I have MM5 model forecasted precip on one curvilinear grid. The NetCDF
> highlights are:
>
> float rain_non(time, i_cross, j_cross) ;
> rain_non:long_name = "ACCUMULATED NONCONVECTIVE PRECIPITATION" ;
> rain_non:units = "cm" ;
> float latitcrs(i_cross, j_cross) ;
> latitcrs:long_name = "LATITUDE (SOUTH NEGATIVE)" ;
> latitcrs:units = "degrees" ;
> float longicrs(i_cross, j_cross) ;
> longicrs:long_name = "LONGITUDE (WEST NEGATIVE)" ;
> longicrs:units = "degrees" ;
>
> And I have NCEP stage4 precip data on a different curvilinear grid. Again,
> the NetCDF highlights.
>
> float APCP_SFC(YCells, XCells) ;
> APCP_SFC:long_name = "Total precipitation" ;
> APCP_SFC:units = "inch" ;
> float longitude(YCells, XCells) ;
> longitude:long_name = "longitude" ;
> longitude:units = "degrees_east" ;
> float latitude(YCells, XCells) ;
> latitude:long_name = "latitude" ;
> latitude:units = "degrees_north" ;
>
>
> So, thinking like a Ferret, I would think I need to regrid these two
> curvilinear grids onto a common regular lat-lon grid and then proceed with
> the comparison between the two. Unfortunately, I can't figure out how to
> do this :( From the docs, I see this in Chp 8 Sec 8.2
>
> "To perform the analysis in a rectilinear coordinate system, the
> conversion of the curvilinear data is most easily done with SAMPLEXY_CURV
> (under development06/00 )."
>
> I can't seem to find that samplexy_curve was implemented. So I searched
> the email archives and found this post which suggests to use scat2grid
> functions.
>
> http://ferret.pmel.noaa.gov/Ferret/Mail_Archives/fu_2001/msg00084.html
>
> But those functions don't seem to work for me. The resulting grid is
> always "No valid data". For example I did:
>
> yes? use APCP_06171200.nc ! My NCEP data
> yes? DEFINE AXIS/X=100w:90w:0.5/UNITS=degrees x0
> yes? DEFINE AXIS/Y=40n:50n:0.5/UNITS=degrees y0
> yes? show data
> currently SET data sets:
> 1> ./APCP_06171200.nc (default)
> name title I J K L
> APCP_SFC Total precipitation 1:1121 1:881 ... ...
> MAPPROJECTION
> MapProjection ... ... ... ...
> LONGITUDE
> longitude 1:1121 1:881 ... ...
> LATITUDE latitude 1:1121 1:881 ... ...
>
> yes? LET outg = SCAT2GRIDGAUSS_XY(longitude, latitude, APCP_SFC, x[gx=x0],
> y[gy=y0],0.5,0.5,2,2)
> yes? shade outg
>
>
> I am using Ferret 5.53 on Fedora Core 1 and am probably doing something
> ignorant :(
>
> thanks,
> daryl
>
--
___________________________________________________
Jaison Kurian
Centre for Atmospheric and Oceanic Sciences
Indian Institute of Science
B A N G A L O R E 560 012
Ph: +91-80-3942505
+91-80-3600450
Fax:+91-80-3600865
___________________________________________________
\ cancel mode verify
!
! curv2regular_demo.jnl
!
! 1. get some sample data on curvilinear co-ordinates
! ---------------------------------------------------
! like
! lat(ni,nj), lon(ni,nj), var(ni,nj,time)
! define ni & nj axis ; suppose we need
! 10 lon points
! 12 lat points &
! 10 time points
! define two "index" axes and a time axis
define axis/x=1:10:1 niax
define axis/y=1:12:1 njax
define axis/t="01-jan-2003:12":"10-jan-2003:12":1/t0="31-dec-2002:12"/np=10/units=days tax
! get curvilinear lat & lon : curv_*
let curv_lon = x[gx=niax] + y[gy=njax]*0
let curv_lat = x[gx=niax]*0 + y[gy=njax]
! get curvilinear variable : curv_var
! "curv_var" is defined in the same way as given in the
! objective_analysis_demo.jnl but with 'time' dimension
let xpts = x[gx=niax] ; let ypts = y[gy=njax]
let wave = sin(kx*xpts + ky*ypts - phase) / 3
let phase = 0 ; let kappa = 0.4
let kx = 0.4 ; let ky = 0.7
let fcn1 = sin(r)/(r+1)
let r = ((xpts-x0)^2+ 5*(ypts-y0)^2)^0.5
let x0 = 3 ; let y0 = 8
let curv_dummy = fcn1 + wave
! to make curv_dummy a 3d variable with time dimension
! like curv_var(ni,nj,time)
let curv_var = curv_dummy + sin(t[gt=tax])
fill/l=10/title="varialbe on curvilinear grid" curv_var
say ; message "let us do the regridding to regular lat-lon grid"
say
! 2. get the 2D lat and lon to 1D ( or to "flatten" lat and lon) : scat_*
!--------------------------------------------------------------------------
let scat_lon = xsequence(curv_lon)
let scat_lat = xsequence(curv_lat)
! 3. Define a regular grid/axes on which the data is needed
!------------------------------------------------------------
! get the safe minimum and maximum for the regular grids to be defined
let x_min = int(scat_lon[i=@min]) ; let x_max = int(scat_lon[i=@max])
let y_min = int(scat_lat[i=@min]) ; let y_max = int(scat_lat[i=@max])
! define axes of regular grid : reg_*
define axis/x=`x_min`:`x_max`:1/units=longitudes reg_lon
define axis/y=`y_min`:`y_max`:1/units=latitudes reg_lat
! 4. Do actual "REGRIDDING" : for details refer Ch3 Sec2.3.35. SCAT2GRIDGAUSS_XY
! (SCAT2GRIDLAPLACE_XY also can be used)
!---------------------------------------------------------------------------------
! define arguments for scat2gridgauss_xy
define symbol scale = 0.5
define symbol cutoff = 2
let reg_var = scat2gridgauss_xy(scat_lon,scat_lat,curv_var,x[gx=reg_lon],y[gy=reg_lat],($scale),($scale),($cutoff),($cutoff))
! 5. Do calculations / plotting or save to a netcdf file
!--------------------------------------------------------
set var/title="Variable on regular grid"/units="myunits" reg_var
fill/l=10 reg_var
!save/file=reg_var.nc/append reg_var
! end of curv2regular_demo.jnl
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement