Hi everyone,
I am trying to plot data that I have on a polar stereographic projection.
The projection info is :
X Position of north pole in grid point: 6.0
Y Position of north pole in grid point: 65.0
Angle between X axis and Greenwich Meridian: 21.0
Projection TRUE at 60 deg N
Grid space at 60 deg N: 100.0 km
I generated 2D latitude and longitude vectors.
Then I used the instructions on how to plot native stereographic plot with the
following code. The problem is that I cannot find the relative center longitude
correctly. Also it seems that if this longitude is not in the region I am
plotting, it doesn't work. Any help? Thanks.
external FNOM "/io/page/ncl/libnclrpn.so"
external FSTOUV "/io/page/ncl/libnclrpn.so"
external FCLOS "/io/page/ncl/libnclrpn.so"
external FSTFRM "/io/page/ncl/libnclrpn.so"
external FSTLIR "/io/page/ncl/libnclrpn.so"
external GET_DIM "/io/page/ncl/libnclrpn.so"
external GET_STEP_NIV "/io/page/ncl/libnclrpn.so"
external GET_TYPGRILLE "/io/page/ncl/libnclrpn.so"
external GET_PROJ "/io/page/ncl/libnclrpn.so"
;*****************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;*****************************************************
begin
;*** open file
iun = 1
ier = FNOM::fnom(iun, "/data/cartes/sea/reg17n.00z", "RND", 0)
ier = FSTOUV::fstouv(iun, "RND")
;*** Get ni and nj (grid dimensions)
ni = 0
nj = 0
GET_DIM::get_dim(iun, "GZ", " ", ni, nj)
;*** Get grid type
grtyps = new(1,character)
GET_TYPGRILLE::get_typgrille(iun, "GZ", " ", grtyps)
grtyp = chartostring(grtyps)
;*** Get latitude and longitude
lon = new((/ni,nj/),float)
lat = new((/ni,nj/),float)
dgrw = 0.0
GET_PROJ::get_proj(iun, "GZ", " ", lon, lat, dgrw, ni, nj)
;* Compute projection center longitude
centerlon = (dgrw-90.0) * 180.0 / 90.0
;************************************************
; create plot
;************************************************
wks = gsn_open_wks("ps","test_map") ; open an ps plot
res = True ; plot mods desired
res@cnFillOn = False
res@cnLinesOn = True
res@gsnAddCyclic = False ; regional data
res@cnLevelSelectionMode = "ManualLevels" ; manual levels
res@cnLevelSpacingF = 4.
res@tiMainString = "Native Stereographic Example"
res@gsnMaximize = True
; the following resources are REQUIRED to plot this projection correctly
res@mpProjection = "Stereographic" ; projection
res@mpLimitMode = "Corners" ; method to zoom
res@mpLeftCornerLatF = lat(0,0)
res@mpLeftCornerLonF = lon(0,0)
res@mpRightCornerLatF = lat(ni-1,nj-1)
res@mpRightCornerLonF = lon(ni-1,nj-1)
res@mpRelativeCenterLon = True ; set a center lon
res@mpCenterLonF = centerlon
res@mpRelativeCenterLat = True ; set a center lat
res@mpCenterLatF = 90. ; center lat
res@tfDoNDCOverlay = True ; do not transform data
;*** Get wanted data
ip1 = 0
ip2 = 12
ip3 = -1
n3 = 1
buffer = new((/ni,nj/),float)
ier = FSTLIR::fstlir(buffer, iun, ni, nj, n3, -1, " ", ip1, ip2, ip3, \
" ", "PN")
buffer@units = "hPa"
buffer@long_name = "mean sea-level pressure"
;*** Plot
plot = gsn_contour_map (wks,buffer,res)
;*** Close file
ier = FSTFRM::fstfrm(iun)
FCLOS::fclos(iun)
end
Christian Page
page AT sca.uqam.ca http://meteocentre.com/
Agent de recherche et de planification +1 514 987 3000 ext. 2376
Departement des Sciences de la Terre et de l'Atmosphere UQAM
This archive was generated by hypermail 2b29 : Tue Feb 19 2002 - 08:48:13 MST