Re: [ferret_users] Projected data

To: ferret_users@xxxxxxxx
Subject: Re: [ferret_users] Projected data
From: "Ingo Bethke (uib)" <Ingo.Bethke@xxxxxxxxxxxxxx>
Date: Tue, 18 Oct 2005 20:58:12 +0800
In-reply-to: <20051017080333.84333.qmail@web32813.mail.mud.yahoo.com>
References: <20051017080333.84333.qmail@web32813.mail.mud.yahoo.com>
Sender: owner-ferret_users@xxxxxxxxxxxxx
User-agent: Mozilla Thunderbird 1.0.7-1.1.fc4 (X11/20050929)

Please find attached an example (I've tried to find some data which is similar to Tonys data). There is still a small but fixable problem where longitudes switch from ~360 to ~0. The sample data can be download from ftp://ftp.gfi.uib.no/ingo/ferret/sample.nc.bz2
Cheers, Ingo

jagadish karmacharya wrote:

Hi,
will you please write and illustrative example.
Thanks
Jagadish
--- "Ingo Bethke (uib)" <Ingo.Bethke@student.uib.no>
wrote:

Hi,

If you are not happy with the default result of the
3-argument versions of the SHADE, FILL , and CONTOUR commands you can
try the POLYGON command in combination with the in ferret available
map projections.

I have tested this on curvlinear model output where
the cell-corners can be used directly to define the polygons. However, it
shouldn't be to difficult to define pseudo-cells (I guess that is
exactly what is done in the SHADE command) which suit your data. The
attached figures crudely illustrate how the result from the curvlinear output
looks like.

Following steps are involved:
The polygon-coordinates need to be gathered in
longitude and a latitude matrices where the first dimension corresponds to
the number of points in each polygon and the second to your data
dimension. It might be neccessary to check and fix the longitude
information of each polygon (e.g. 179 and -179 -> 179 and 181 ect.). After
calling the projection script (e.g. go mp_orthographic 0 90) you need to
set the variabels mp_x and mp_y to your polygon longitude and latitude
coordinates. Prior to the call of the 3-argument version of the polygon
command (e.g. polygon x_page*mp_mask,y_page*mp_mask,sst) you might also
want to call mp_aspect.

The drawback of this method is that not all
projection scripts are robust to user-defined mp_x,mp_y variables. With the
exception of mp_orthographic I'm still trying hard to get grid
lines and land masks right for the different projections.

Hope this helps...

Cheers, Ingo

Ansley Manke wrote:

Hi Tony,
You can plot this data with the curvilinear forms

of the plotting

commands. The variables LONGITUDE and LATITUDE

are 2-dimensional

curvilinear coordinate fields, so you can use this

command to get the

correct longitude and latitude labellings

fill temperature, longitude, latitude

The fact that the SHOW GRID has the axes of the

variables like

temperature labeled as longitude and latitude is

misleading; the axes

are just indices 1 to 792 and the coordinates are

in the curvilinear

coordinate variables. Look at the SHOW GRID

output:

yes? show grid temperature
GRID GPI1
name axis # pts start

end

X LONGITUDE 792 r 1E

72E(792)

Y LATITUDE 792 r 1N

792N

...

It has a latitude 1 N to 792N, and longitude 1 to

792 as well, which

is clearly incorrect. The curvilinear graphics

commands require that

the X and Y axes of the variable field such as

TEMPERATURE or

SALINITY, and the X and Y axes of the coordinates

LONGITUDE and

LATITUDE agree and then it can draw the plot.

Ansley

Tony Jolibois wrote:

Hi all,

I have to plot some stereographic projected data

(latitude 55N to

90N, projected in the netCDF), and I don't know

how I can use it with

longitude and latitude (instead of x and y).

Here is the ferret command I use to produce the

plot attached. The

axes are not good and I can not fill the land

mask...

[console@rdp1-jaune arc]$ ferret -nojnl
NOAA/PMEL TMAP
FERRET v5.50
Linux 2.4.3-12smp - 01/15/03
12-Oct-05 09:40

yes? use

mercatorPsy3v1R1v_arc_mean_20050923_R20051005.nc

yes? show data
currently SET data sets:
1>

./mercatorPsy3v1R1v_arc_mean_20050923_R20051005.nc (default)

name title I

J

K L
TEMPERATURE
temperature 1:792

1:792

1:43 ...
SALINITY salinity 1:792

1:792

1:43 ...
U zonal velocity 1:792

1:792

1:43 ...
V meridian velocity 1:792

1:792

1:43 ...
KZ vertical eddy diffusivity 1:792

1:792

1:43 ...
LONGITUDE
longitude 1:792

1:792

... ...
LATITUDE latitude 1:792

1:792

... ...
SSH sea surface height 1:792

1:792

... ...
MLD temperature ocean mixed layer t 1:792

1:792

... ...
MLP density ocean mixed layer thick 1:792

1:792

... ...
TAUX windstress eastward Tx componen 1:792

1:792

... ...
TAUY windstress northward Ty compone 1:792

1:792

... ...
QTOT total net heat flux 1:792

1:792

... ...
EMP water flux 1:792

1:792

... ...
QSR surface downward solar heat flu 1:792

1:792

... ...

yes? show grid temperature
GRID GPI1
name axis # pts start

end

X LONGITUDE 792 r 1E

72E(792)

Y LATITUDE 792 r 1N

792N

DEPTH DEPTH (m) 43 i- 0

5500

normal T
yes? set region/x=1:792/y=1:792/z=0
yes? fill/levels=(-3,15,0.2) temperature

Best regards,
Tony

--
Tony Jolibois

CLS - Direction Océanographie Spatiale
8-10 rue hermès - 31520 Ramonville Saint-Agne
tjolibois@cls.fr

=== message truncated ===

__________________________________ Yahoo! Music Unlimited Access over 1 million songs. Try it free.
http://music.yahoo.com/unlimited/

! example script for the use of map projections with curvlinear data 
! (2005,ingo@nersc.no) 


! access data
use sample.nc


! compute cell corners 
let lon_corner_1=0.25*ysequence(lon[i=@shf:-1,j=@shf:-1]\
                 +lon[i=@shf:-1]+lon[j=@shf:-1]+lon)
let lon_corner_2=0.25*ysequence(lon[i=@shf:+1,j=@shf:-1]\
                 +lon[i=@shf:+1]+lon[j=@shf:-1]+lon)
let lon_corner_3=0.25*ysequence(lon[i=@shf:+1,j=@shf:+1]\ 
                 +lon[i=@shf:+1]+lon[j=@shf:+1]+lon)
let lon_corner_4=0.25*ysequence(lon[i=@shf:-1,j=@shf:+1]\
                 +lon[i=@shf:-1]+lon[j=@shf:+1]+lon)

let lat_corner_1=0.25*ysequence(lat[i=@shf:-1,j=@shf:-1]\
                 +lat[i=@shf:-1]+lat[j=@shf:-1]+lat)
let lat_corner_2=0.25*ysequence(lat[i=@shf:+1,j=@shf:-1]\
                 +lat[i=@shf:+1]+lat[j=@shf:-1]+lat)
let lat_corner_3=0.25*ysequence(lat[i=@shf:+1,j=@shf:+1]\ 
                 +lat[i=@shf:+1]+lat[j=@shf:+1]+lat)
let lat_corner_4=0.25*ysequence(lat[i=@shf:-1,j=@shf:+1]\
                 +lat[i=@shf:-1]+lat[j=@shf:+1]+lat)


! concatenate cell corners
let lon_corner_12=xcat(lon_corner_1,lon_corner_2)
let lon_corner_123=xcat(lon_corner_12,lon_corner_3)
let lon_corners=xcat(lon_corner_123,lon_corner_4)

let lat_corner_12=xcat(lat_corner_1,lat_corner_2)
let lat_corner_123=xcat(lat_corner_12,lat_corner_3)
let lat_corners=xcat(lat_corner_123,lat_corner_4)


! mask bad corners (ideally these should be replaced...) 
let mask_corners=if (abs(ysequence(lon)-lon_corners) lt 10) then 1 


! select a time slice and reshape data
set variable /bad=0 fice
let field=ysequence(fice[l=1])


! run projection script 
go mp_orthographic -20 70 


! draw global basemap 
use levitus_climatology 
let dummy=x[g=temp]+y[g=temp]
set grid temp
go mp_aspect
fill /nolab/noaxis/nokey/palette=blue dummy,x_page,y_page
go mp_fland 20 greyscale overlay detailed 
cancel data levitus_climatology.cdf


! plot polygons 
let mp_x=lon_corners*mask_corners
let mp_y=lat_corners*mask_corners
poly /ov/nolab/nokey/noaxis/palette=bluescale/levels=20 \
     x_page*mp_mask,y_page*mp_mask,field


! draw graticules
use levitus_climatology 
set grid temp
go mp_graticule
cancel data levitus_climatology.cdf


! print result
!frame /file=sample.gif

Previous by thread: Re: [ferret_users] Projected data
Next by thread: [ferret_users] Re: ferret -script and arguments

[Thread Prev][Thread Next][Index]

Dept of Commerce / NOAA / OAR / PMEL / TMAP