Hi,
The third argument to the CURV_TO_RECT_MAP function is a 2-dimensional variable, defined on the output grid that you want.
You could use
yes? define axis/x=-5.67:5.25:0.125/units=degrees_east xax
yes? define axis/y=-2:2/units=degrees_north yax
yes? let lonlatout = x[gx=xax] + y[gy=yax]
That last command defines a variable on the rectilinear grid of the result. Its values aren't ever used; the function is just getting that variable's grid in the x and y directions.
On 10/8/2010 10:22 AM, Szymon Roziewski wrote:I need to use (I think I need to) the function CURV_TO_RECT_MAP and CURV_TO_RECT (I want to transform wind data from one bigger curvulinear grid to another smaller rectilinear and next use it in a model)
there is a tutorial on the internet:
yes? ! For convenience define variables with the input grid
yes? let lonin = geolon[d=1]
yes? let latin = geolat[d=1]
yes? ! Define output grid and a variable on the output grid
yes? define axis/x=0:360:5/modulo/units=degrees xax
yes? ! Compute the mapping to the rectangular grid, save to a file
yes? def axis/y=0:85:5/units=degrees yax
yes? let lonlatout = y[gy=yax] + x[gx=xax]
yes? let map = curv_to_rect_map ( lonin,latin,lonlatout,10)
yes? save/clobber/file=curv_map.nc map
so what I need it's a properly lonlatout matrix which I would like to cover the following area:
(bottom left) -5.67E, -2S (bottom right) 5.25W, -2S
(top left) -5.67E, 10N (top right) 5.25W, 10N
So I need one matrix and I don't excatly know how to get it. e.g. in matlab one may get matrices X and Y but curv_to_rect_map demand only one.
2010/10/8 Karl Smith <Karl.Smith@xxxxxxxx>
Hi Simon,
This defined a grid which can be used as the grid for other data. If you want the equivalent of meshgrid as you mentioned in your later e-mail,
set grid grid_out
def xvals = X + 0.0 * Y
def yvals = 0.0 * X + Y
The xvals and yvals are then variables that can be saved or used as desired.
-- Karl
On 10/8/2010 9:15 AM, Szymon Roziewski wrote:
Hi Karl,2010/10/8 Karl Smith <Karl.Smith@xxxxxxxx <mailto:Karl.Smith@xxxxxxxx>>
Thank You for your response,
I tried as You told me and what happened next:
**ERROR: variable unknown or not in data set: GRID_OUT
LIST/FORMAT=CDF/format=(448F8.2)/file=GRIDOUT/clobber grid_out
Command file, command group, or REPEAT execution aborted
STOP -script mode, ERROR RUNNING SCRIPT
ferret doesn't know about values on this grid maybe because treats this
as an abstract object.
Hi Simon,
Using the E/W on the X axis and N/S on the Y axis will identify them
as longitude/latitude. So you were close.
-- Karl
define axis /x=5.666666667E:5.25W:0.083333333 xax
def axis /y=2.0S:10.0N:0.083333333 yax
def grid /x=xax /y=yax grid_out
On 10/8/2010 7:54 AM, Szymon Roziewski wrote:
Now I'm thinking how to make appropriate grid for that:
My data are on grid;
(bottom left) -5.67E, -2S (bottom right) 5.25W, -2S
(top left) -5.67E, 10N (top right) 5.25W, 10N
Maybe there is special function wchich dealing with it in ferret?
2010/10/8 Szymon Roziewski <szymon.roziewski@xxxxxxxxx
<mailto:szymon.roziewski@xxxxxxxxx>
<mailto:szymon.roziewski@xxxxxxxxx
<mailto:szymon.roziewski@xxxxxxxxx>>>
Ok,
I found it,
the problem was about statement:
let grid_out = y[gy=yax] + x[gx=xax]
which I didn't get well.
Cheers
2010/10/8 Szymon Roziewski <szymon.roziewski@xxxxxxxxx
<mailto:szymon.roziewski@xxxxxxxxx>
<mailto:szymon.roziewski@xxxxxxxxx
<mailto:szymon.roziewski@xxxxxxxxx>>>
I also attach my script which does compute coefficient
matrix.
DEFINE AXIS/x=1:448:1/unit=degree xaxis
DEFINE AXIS/y=1:615:1/unit=degree yaxis
DEFINE GRID/x=xaxis/y=yaxis gridlonlat
!read curvilinear coordinates onto grid 448x615
FILE/VARIABLES=clon_in/COLUMNS=448/GRID=gridlonlat
"/home/szymon/WAM/CHECK/umlonicm448x615.dat"
FILE/VARIABLES=clat_in/COLUMNS=448/GRID=gridlonlat
"/home/szymon/WAM/CHECK/umlaticm448x615.dat"
!create grid which we desire
define
axis/x=-5.666666667:5.25:0.083333333/modulo/units=degrees xax
def axis/y=-2.0:10.0:0.083333333/units=degrees yax
let grid_out = y[gy=yax] + x[gx=xax]
!do transformation matrix
let map = CURV_TO_RECT_MAP(clon_in[d=1], clat_in[d=2],
grid_out, 4)
!save data
save/clobber/file=curv_map_UM_TO_ICM.nc map
and another script which uses it:
DEFINE AXIS/x=1:448:1/unit=degree xaxis
DEFINE AXIS/y=1:615:1/unit=degree yaxis
DEFINE GRID/x=xaxis/y=yaxis gridlonlat
!load data for transformation - wind fields u,v
file/var=wlon/grid=gridlonlat/format=stream/type=r4/swap
"/home/szymon/WAM/CHECK/03225_2010092100+48000000c3hs000000000000000"
file/var=wlat/grid=gridlonlat/format=stream/type=r4/swap
"/home/szymon/WAM/CHECK/03226_2010092100+48000000c3hs000000000000000"
!using transform matrix
use "/home/szymon/ferret/vectors/curv_map_UM_TO_ICM.nc"
! do transformation
let wlon_icm = curv_to_rect(wlon[d=1], map[d=3])
let wlat_icm = curv_to_rect(wlat[d=2], map[d=3])
!save transformation data which You can see for example
in the
picture
*wind_2010092103_ICM_remapped_132x145_and_interpolated.gif* or
*UM_Significant_Wave_Height_21_September_2010_020000_TMP.gif
(vectors only)*
save/format=(448F8.2)/file=winduumicm+48.txt/clobber
wlon_icm
save/format=(448F8.2)/file=windvumicm+48.txt/clobber
wlat_icm
Best regards,
Simon
2010/10/8 Szymon Roziewski <szymon.roziewski@xxxxxxxxx
<mailto:szymon.roziewski@xxxxxxxxx>
<mailto:szymon.roziewski@xxxxxxxxx
<mailto:szymon.roziewski@xxxxxxxxx>>>
Hello there,
I have problem with appropriate remapping from
curvilinear
coordinates to rectilinear coordinates.
The area which is covered with curvilinear
coordinates is
bigger than I need in rectilinear coordinates (it's
a subset
of curvilinear grid).
I would like to raise such a question if ferret can
do this
transformation automaticly i.e. I needn't to care
about the
curvilinear grid is about bigger area than rectilinear.
Or I need to tell ferret that fact by some
computations and
findings.
I attach files:
- wind_2010092103_ICM_aspect.gif - this is wind on
curvilinear grid it is a bigger area
-
wind_2010092103_ICM_remaped_132x145_and_interpolated.gif -
this is wind on rectilinear grid it's a smaller area.
As You can see ferret did transformation but in a way I
didn't intend. The Baltic Sea area is a little bit
translated to the north and the northern part of sea was
cut. The picture should wrap all area of Baltic Sea.
In the picture
UM_Significant_Wave_Height_21_September_2010_020000_TMP.gif
I marked the area which matches interpolated wind
field on
Baltic Sea and You can see that it disagrees to
significant
wave height which occupies properly area. It means that
interpolation from curvilinear grid 448x615 to the
rectilinear grid 132x145 was done improperly.
I hope I was understood.
Kind regards,
Szymon Roziewski
--
Karl M. Smith, Ph.D.
Research Scientist and Software Developer
UW/JISAO and NOAA/PMEL/TMAP
karl.smith@xxxxxxxx
(206) 526-4806
"The contents of this message are mine personally and do
not necessarily reflect any position of the Government
or the National Oceanic and Atmospheric Administration."
--
Z wyrazami szacunku,
Szymon Roziewski