Plotting REMO data in ncl

Dear all,

with the hint of Axel I could plot REMO data on a Cylindrical
Equidistant grid.

In case you want to do the same here is what I did ;)

1.) derotating the rlat / rlon values out of the netcdf file with the
rotated north pole (see code section 2)
2.) setting the map center to the rotated south pole (with the
transformation of of Axel)

Cheers,
Christof

CODE:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

;------------------------------------
; 1. data & matadata
;------------------------------------
; the file to open
;------------------------------------
file1 = addfile ("t1_167_e040009m1989_2008.nc", "r")

;------------------------------------
; getting lat&lon information
;------------------------------------
rotlat= file1->rlat
rotlon= file1->rlon

;----------------------------------------
; getting testdata to plot something
;----------------------------------------
P = file1->var167

dims_P = getfilevardims(file1,"var167")

p = P($dims_P(0)$|:,$dims_P(1)$|:,$dims_P(2)$|:,$dims_P(3)$|:)

;--------------------------------------
; getting dimensions of nlat and nlon
;--------------------------------------
nlat = dimsizes(rotlat)
nlon = dimsizes(rotlon)

;----------------------------------------
; coordinates of the rotated north pole
;----------------------------------------

n_pollat = 42.5
n_pollon = 83.0

;----------------------------------------
; coordinates of the rotated south pole
;----------------------------------------

s_pollat = -42.5
s_pollon = -97.0

;---------------------------------------------------
; 2. derotation of the lat/lon information
;---------------------------------------------------
;
; to derotate the lat/long information
; for setting the rigth "corners" in res@mpLimitMode
;
;----------------------------------------------------

derot_ll_lon = new( (/ nlon, nlat /),double)
derot_ll_lat = new( (/ nlat, nlon /),double)

;----------------------------------------------------
; getting the rotated longitudes & latitudes
;----------------------------------------------------
pi = 4.*atan(1.)
pi2deg = 180./pi
deg2pi = pi/180

rsinpol = sin(deg2pi*n_pollat)
rcospol = cos(deg2pi*n_pollat)

;----------------------------------------------------
; derotate the longitudes
;----------------------------------------------------

do i=0, nlon-1
  do j=0, nlat-1
 
    derot_ll_lon(i,j) = rotlon(i)

    if( rotlon(i) .GT. 180.) then
      derot_ll_lon(i,j) = rotlon(i) - 360.
    end if

    if ( rotlon(i) .LT. -180. ) then
      derot_ll_lon(i,j) = rotlon(i) + 360.
    end if

    help_rotlon = deg2pi * derot_ll_lon(i,j)
    help_rotlat = deg2pi * rotlat(j)

    rarg1 = -sin(help_rotlon)
    rarg2 = rcospol * tan(help_rotlat) - rsinpol * cos(help_rotlon)

    if ( rarg2 .EQ. 0 ) then
      rarg2 = 1.E-20
    end if

    derot_ll_lon(i,j) = n_pollon + pi2deg * atan2(rarg1,rarg2)

    if ( derot_ll_lon(i,j) .LT. -180. ) then
      derot_ll_lon(i,j) = derot_ll_lon(i,j) + 360
    end if

    if ( derot_ll_lon(i,j) .GT. 180. ) then
      derot_ll_lon(i,j) = derot_ll_lon(i,j) - 360
    end if

  end do
end do

; --------------------------------------
; derotate the latitudes
; --------------------------------------

    pi = 4.*atan(1.)
    pi2deg = 180./pi
    deg2pi = pi/180

    rsinpol = sin(deg2pi*n_pollat)
    rcospol = cos(deg2pi*n_pollat)

do i=0, nlat-1
  do j=0, nlon-1

    help_lon = rotlon(j)

    if(rotlon(j) .GT. 180.) then
      help_lon = rotlon(j) - 360.
    end if

    if(rotlon(j) .LT. -180.) then
      help_lon = rotlon(j) + 360.
    end if

    help_lon = deg2pi * help_lon
    help_lat = deg2pi * rotlat(i)

    rarg = rsinpol * sin(help_lat) + rcospol * cos(help_lon) * cos(help_lat)

    derot_ll_lat(i,j) = pi2deg * asin(rarg)

  end do
end do

;----------------------------------------------------
; 3. set the workspace
;----------------------------------------------------

wks = gsn_open_wks ("ps", "lcnative")
gsn_define_colormap (wks,"gui_default")

;---------------------------------------
; 4. set up map & plot
;---------------------------------------
res = True ; plot mods desired

res@mpProjection = "CylindricalEquidistant" ; projection

;-------------------------------
; configuring the map center
;-------------------------------
res@mpCenterLonF = s_pollon
res@mpCenterLatF = s_pollat+90.0

res@tfDoNDCOverlay = True

res@cnFillOn = True ; color fill
res@cnLinesOn = False ; no contour lines
res@cnLineLabelsOn = False ; no contour labels
res@gsnSpreadColors = True ; use total colormap
res@gsnSpreadColorStart = 4
res@gsnSpreadColorEnd = -1
res@cnInfoLabelOn = False ; no contour info label

res@tmXBOn = False ; tickmarks are wrong anyway
res@tmXTOn = False
res@tmYLOn = False
res@tmYROn = False

res@mpGridLineDashPattern = 2 ; lat/lon lines as dashed
res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks
 
;----------------------------
; Title
;----------------------------
res@tiMainString = "The method to plot REMO data"
res@tiMainFontHeightF = 0.020 ; smaller title

;----------------------------
; Corners
;----------------------------
res@mpLimitMode = "Corners" ; choose range of map
res@mpLeftCornerLatF = derot_ll_lat(0,0)
res@mpLeftCornerLonF = derot_ll_lon(0,0)
res@mpRightCornerLatF = derot_ll_lat(nlat-1,nlon-1)
res@mpRightCornerLonF = derot_ll_lon(nlon-1,nlat-1)

;---------------------------
; Plot
;---------------------------
plot = gsn_csm_contour_map(wks,p(0,0,:,:),res)

end

Axel Seifert wrote:
> Dear Christof,
>
> in NCL a rotated lat/lon is a "CylindricalEquidistant" projection. The
> two code snippets below should help you to get either the native
> projection or transform to another projection. In the latter case you
> need the lat/lon in 2d arrays. The first method gets the grid
> information (corner points and pole) from a GRIB file, but NetCDF should
> provide the same information.
>
> ; for native grid plot
>
> res@mpProjection = "CylindricalEquidistant" ; projection
> res@mpLimitMode = "Corners" ; method to zoom
> res@mpLeftCornerLatF = b->g10_lat_0@corners(0)
> res@mpLeftCornerLonF = b->g10_lon_1@corners(0)
> res@mpRightCornerLatF = b->g10_lat_0@corners(2)
> res@mpRightCornerLonF = b->g10_lon_1@corners(2)
> res@mpCenterLonF = b->g10_lon_1@Longitude_of_southern_pole
> res@mpCenterLatF = b->g10_lon_1@Latitude_of_southern_pole + 90.
> res@tfDoNDCOverlay = True ; dont transform data
>
> res@tmXBOn = False ; tickmarks are wrong anyway
> res@tmXTOn = False
> res@tmYLOn = False
> res@tmYROn = False
>
> ; for regular lat/lon projection
>
> resLatLon@trGridType = "Curvilinear"
> resLatLon@sfXArray = lon2d
> resLatLon@sfYArray = lat2d
> resLatLon@mpLimitMode = "LatLon"
> resLatLon@tfDoNDCOverlay = False
> resLatLon@pmLabelBarOrthogonalPosF = 0.20
>
> Hope that helps,
> Axel
>
>
>
>
>
> Am Mittwoch, den 20.01.2010, 14:52 +0100 schrieb Christof Wilhelm:
>
>> Dear ncl-talk readers,
>>
>> is there a possibility in ncl to plot geographical coordinates in a
>> spherical coordinate system?
>>
>> The background is the following, that I generally want to find a method
>> to plot REMO output with ncl.
>> REMO output comes on rotated geographical coordinates.
>>
>> I tried it with the method of native grids following the link below, but
>> couldn't solve the problem, because the REMO data are not projected in
>> any of the common ncl projections.
>>
>> http://www.ncl.ucar.edu/Support/talk_archives/2009/0052.html
>>
>>
>> THX for any help,
>> Christof
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
>

-- 
Christof Wilhelm 
Max-Planck-Institute for Meteorology
Bundesstraße 53
20146 Hamburg
Germany
T +49 40 41173 433
M christof.wilhelm@zmaw.de
---------------------------------------------------------------------------
http://www.mpimet.mpg.de/en/institut/mitarbeiter/wilhelmchristof/index.html
---------------------------------------------------------------------------
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk