Personal tools

##### Sektionen
You are here: Home / Rotated Grid
Info
Our docs have moved to https://docs.dkrz.de This portal will be shutdown as soon the last content is migrated.

# Rotated Grid

How to plot data computed on a rotated grid on the unrotated grid location. Here, CORDEX EUR-11 data are used to demonstrate the transform functions inside the example script.

Example data file (CORDEX):  tas_rotated_grid_EUR11.nc

Example script:

```;------------------------------------------------------------
;-- set global constants
;------------------------------------------------------------
pi       =  4.0*atan(1.)
fillval  = -99999.9

;------------------------------------------------------------
;-- Function:      create_lon2d(y,x)
;-- Description:   create 2d lon array from 1D rotlon, rotlat
;------------------------------------------------------------
undef("create_lon2d")
function create_lon2d(rotlat[*]:numeric, rotlon[*]:numeric)
local x, i
begin
x = new((/dimsizes(rotlat),dimsizes(rotlon)/),typeof(rotlat))
do i=0,dimsizes(rotlat)-1
x(i,:) = rotlon
end do
return(x)
end

;------------------------------------------------------------
;-- Function:      create_lat2d(y,x)
;-- Description:   create 2d lat array from 1D rotlon,rotlat
;------------------------------------------------------------
undef("create_lat2d")
function create_lat2d(rotlat[*]:numeric, rotlon[*]:numeric)
local y, i
begin
y = new((/dimsizes(rotlat),dimsizes(rotlon)/),typeof(rotlat))
do i=0,dimsizes(rotlon)-1
y(:,i) = rotlat
end do
return(y)
end

;------------------------------------------------------------
;-- Function:      unrot_lon(rotlat,rotlon,pollat,pollon)
;-- Description:   transform rotated longitude to longitude
;------------------------------------------------------------
undef("unrot_lon")
function unrot_lon( rotlat:numeric, rotlon:numeric, pollat[1]:numeric, pollon[1]:numeric )
local rotlat, rotlon, pollon, pollat, lon, s1, c1, s2, c2, rlo, rla, i, tmp1, tmp2
begin
lon = fillval
lon@_FillValue = fillval

if (any(dimsizes(rotlat) .ne. dimsizes(rotlon)) .and. \
(dimsizes(dimsizes(rotlat)).ne.1 .or. dimsizes(dimsizes(rotlon)).ne.1)) then
print("Function unrot_lon: unrot_lon:  rotlat and rotlon dimensions do not match")
return(lon)
end if

if (dimsizes(dimsizes(rotlat)).eq.1 .and. dimsizes(dimsizes(rotlon)).eq.1) then
rla = create_lat2d(rotlat,rotlon)    ;-- create 2D latitude array
rlo = create_lon2d(rotlat,rotlon)    ;-- create 2D longitude array
else
rla = rotlat
rlo = rotlon
end if

lon := (/rlo/)                         ;-- reassign lon
lon@_FillValue=fillval

tmp1 = s2*(-s1*cos(rlo)*cos(rla)+c1*sin(rla))-c2*sin(rlo)*cos(rla)
tmp2 = c2*(-s1*cos(rlo)*cos(rla)+c1*sin(rla))+s2*sin(rlo)*cos(rla)

lon@units = "degrees_east"
print("Function unrot_lon: min/max     "+sprintf("%8.4f", min(lon(0,:)))+"  "+sprintf("%8.4f", max(lon(0,:))))

delete([/rlo,rlo,c1,s1,c2,s2,tmp1,tmp2/])

return(lon)
end

;------------------------------------------------------------
;-- Function:      unrot_lat(rotlat,rotlon,pollat,pollon)
;-- Description:   transform rotated latitude to latitude
;------------------------------------------------------------
undef("unrot_lat")
function unrot_lat( rotlat:numeric, rotlon:numeric, pollat[1]:numeric, pollon[1]:numeric )
local rotlat, rotlon, pollon, pollat, lat, s1, c1, rlo, rla, i
begin
lat = fillval
lat@_FillValue = fillval

if (any(dimsizes(rotlat) .ne. dimsizes(rotlon)) .and. \
(dimsizes(dimsizes(rotlat)).ne.1 .or. dimsizes(dimsizes(rotlon)).ne.1)) then
print("Function unrot_lat:  rotlat and rotlon dimensions do not match")
return(lat)
end if

if (dimsizes(dimsizes(rotlat)).eq.1 .and. dimsizes(dimsizes(rotlon)).eq.1) then
rla = create_lat2d(rotlat,rotlon)    ;-- create 2D latitude array
rlo = create_lon2d(rotlat,rotlon)    ;-- create 2D longitude array
else
rla = rotlat
rlo = rotlon
end if

lat := (/rla/)                         ;-- reassign lat
lat@_FillValue=fillval

lat = s1*sin(rla)+c1*cos(rla)*cos(rlo)

lat@units = "degrees_north"
print("Function unrot_lat: min/max     "+sprintf("%8.4f", min(lat(:,0)))+"  "+sprintf("%8.4f", max(lat(:,0))))

delete([/rlo,rla,c1,s1/])

return(lat)
end

;----------------
;--  MAIN
;----------------
begin
;-- open file and read variables
var     =  f->tas
rlat    =  f->rlat
rlon    =  f->rlon
rotpole =  f->rotated_pole
pollat  =  rotpole@grid_north_pole_latitude
pollon  =  rotpole@grid_north_pole_longitude

;-- unrotate the grid and set 2D lat/lons
lon2d      =  unrot_lon(rlat, rlon, pollat, pollon)
lat2d      =  unrot_lat(rlat, rlon, pollat, pollon)
var@lat2d  =  lat2d                                ;-- 2D latitudes
var@lon2d  =  lon2d                                ;-- 2D longitudes

;-- calculate the min and max lat/lons for the map plot
minlat  =  min(lat2d)                              ;-- retrieve minimum latitude value
minlon  =  min(lon2d)                              ;-- retrieve maximum latitude value
maxlat  =  max(lat2d)                              ;-- retrieve minimum longitude value
maxlon  =  max(lon2d)                              ;-- retrieve maximum longitude value

;-- open a workstation
wks = gsn_open_wks("png","plot_rotated_grid")

;-- set resources
res                       =  True
res@gsnFrame              =  False                 ;-- don't advance frame
res@gsnAddCyclic          =  False                 ;-- data are not global, don't add lon cyclic point

res@pmTickMarkDisplayMode = "Always"               ;-- draw nicer tickmarks

res@mpDataBaseVersion     = "MediumRes"            ;-- choose map database
res@mpMinLatF             =  minlat - 1.           ;-- set min lat
res@mpMaxLatF             =  maxlat + 1.           ;-- set max lat
res@mpMinLonF             =  minlon - 1.           ;-- set min lon
res@mpMaxLonF             =  maxlon + 1.           ;-- set max lon
res@mpGridAndLimbOn       =  True                  ;-- turn on grid lines

res@cnFillOn              =  True                  ;-- turn on contour fill
res@cnLinesOn             =  False                 ;-- don't draw contour lines
res@cnFillPalette         = "BlueYellowRed"        ;-- choose color map

res@lbLabelBarOn          =  True                  ;-- turn on labelbar

res@tiMainString          = "NCL Doc: rotated grid"  ;-- title
res@tiMainOffsetYF        = -0.025                 ;-- move title downward

res@vpWidthF              =  0.6                   ;-- width of viewport
res@vpHeightF             =  0.48                  ;-- height of viewport

;-- create the first plot
res@vpXF                  =  0.12                  ;-- start x-position
res@vpYF                  =  1.02                  ;-- start y-Position

plot1 = gsn_csm_contour_map(wks,var(0,0,:,:),res)  ;-- use default projection (CE)

;-- create the second plot
delete(res@tiMainString)                           ;-- we don't need the title twice

res@vpXF                  =  0.15                  ;-- start x-position
res@vpYF                  =  0.493                 ;-- start y-position

res@mpProjection          = "Orthographic"         ;-- change projection
res@mpCenterLatF          =  minlat + (maxlat -minlat)/2 ;-- center point of view latitude
res@mpCenterLonF          =  minlon + (maxlon -minlon)/2 ;-- center point of view longitude
res@mpLimitMode           = "LatLon"               ;-- map limits mode
res@mpMinLatF             =  minlat - 1.           ;-- set min lat
res@mpMaxLatF             =  maxlat + 1.           ;-- set max lat
res@mpMinLonF             =  minlon - 1.           ;-- set min lon
res@mpMaxLonF             =  maxlon + 1.           ;-- set max lon
res@mpPerimOn             =  False                 ;-- don't draw the box around the plot

res@lbOrientation         = "vertical"             ;-- vertical label bar
res@lbLabelStride         =  2
res@lbLabelPosition       = "Left"                 ;-- labelbar labels on left side
res@pmLabelBarOrthogonalPosF = -1.37               ;-- labelbar on the left side

res@tmXTLabelDeltaF       = -0.5                   ;-- decrease space between ticks and labels
res@tmXBLabelDeltaF       = -0.5                   ;-- decrease space between ticks and labels
res@tmYLLabelDeltaF       = -0.5                   ;-- decrease space between ticks and labels
res@tmYRLabelDeltaF       = -0.5                   ;-- decrease space between ticks and labels

plot2 = gsn_csm_contour_map(wks,var(0,0,:,:),res)  ;-- draw second plot

;-- draw text
txres                =  True
txres@txFontHeightF  =  0.016
txres@txJust         = "CenterLeft"

gsn_text_ndc(wks,"Projection:",             0.74, 0.91, txres)  ;-- next to first plot
gsn_text_ndc(wks,"Cylindrical Equidistant", 0.74, 0.89, txres)

gsn_text_ndc(wks,"Projection:",  0.77, 0.43, txres)   ;-- next to second plot
gsn_text_ndc(wks,"Orthographic", 0.77, 0.41, txres)