You are here: Home / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL example rotated grid: mask data by shapefile content (country borderline)

DKRZ NCL example rotated grid: mask data by shapefile content (country borderline)

This example shows how to mask data on a rotated curvilinear grid by a shapefile of Argentina's borderline (download shapefile from It uses the function shapefile_mask_data from shapefile_utils.ncl ( to create a mask array which is used with the mask function. A small trick brings shapefile_mask_data to accept the named 2D coordinates of the rotated grid.

Example script:

   DKRZ NCL example:

           Read 2D-coordinates on a rotated grid and use a shapefile to mask the data
           for a selected region.

           Download the shapefile from
           gadm36_ARG_shp/gadm36_ARG_0.shp     Argentina borderline

   Based on NCL example:
        shapefile_utils.ncl contains functions for shapefiles and can be copied
         The function shapefile_mask_data() is used to mask the data by a shapefile
         containing the borderline of Argentina.
   Known bug:
           For rotated grids or named coordinates of curvilinear grids the function
           shapefile_mask_data() is running into an error because the named
           coordinates let it go into the part of rectilinear grids (bug). To prevent
           running into this error you simply have to multiply the variable by 1. and
           the names of the coordinates will be stripped off. ;-)


;-- load shapefile functions library file
load "shapefile_utils.ncl"

;-- open data file and read the variable
f   = addfile("","r")
var = f->T_CL(0,:,:) * 1.     ;-- read variable and strip off the names of the coordinates
lat = f->lat
lon = f->lon
var@lat2d = lat               ;-- tell NCL to use 2D coordinates
var@lon2d = lon               ;-- tell NCL to use 2D coordinates

;-- shapefile name; downloaded from
shpname = "gadm36_ARG_shp/gadm36_ARG_0.shp"
shpvar  = "NAME_0"

shpf = addfile(shpname,"r")                   ;-- open the shapefile

;-- set sub-region
minlat = min(shpf->y) - 2.5
maxlat = max(shpf->y) + 1.5
minlon = min(shpf->x) - 2.5
maxlon = max(shpf->x) + 1.5

;-- graphics
wks = gsn_open_wks("png","plot_REMO_Argentina_mask_data_by_shapefile_0")

;-- create mask array
opt                 =  True                   ;-- shapefile mask resources
opt@return_mask     =  True                   ;-- this forces the return of a 0s and 1s mask array
opt@debug           =  True                   ;-- show debugging information

mask_array          =  shapefile_mask_data(var,shpname,opt)
mask_array!0        = "y"                     ;-- add name coordinate lat
mask_array!1        = "x"                     ;-- add named coordinate lon
mask_array@lat2d    =  lat                    ;-- commit lat values
mask_array@lon2d    =  lon                    ;-- commit lon values
mask_array@coordinates = "lat lon"

;-- mask the data with the mask_array
mask_data           =  mask(var,mask_array,1) ;-- mask var by mask_array
mask_data!0         = "y"                     ;-- add name coordinate y
mask_data!1         = "x"                     ;-- add named coordinate x
mask_data@lat2d     =  lat                    ;-- commit lat values
mask_data@lon2d     =  lon                    ;-- commit lon values
mask_data@coordinates = "lat lon"

;-- masked contour plot
cnres               =  True
cnres@gsnDraw       =  False
cnres@gsnFrame      =  False

cnres@vpXF          =  0.1
cnres@vpYF          =  0.935
cnres@vpWidthF      =  0.8
cnres@vpHeightF     =  0.8

cnres@cnFillOn      =  True
cnres@cnFillPalette = "BlueYellowRed"

cnres@mpOutlineOn   =  False                  ;-- don't use NCL's map outlines
cnres@mpMinLatF     =  minlat
cnres@mpMaxLatF     =  maxlat
cnres@mpMinLonF     =  minlon
cnres@mpMaxLonF     =  maxlon

cnres@lbBoxMinorExtentF = 0.1                ;-- decrease the labelbar box height
cnres@pmLabelBarOrthogonalPosF = -0.03       ;-- move the labelbar upward

cnres@tiMainString  = "Argentina"

plot = gsn_csm_contour_map(wks,mask_data,cnres) ;-- create contour plot

;-- add shapefile polyline content
plres               =  True
plres@gsLineColor   = "black"
pid = gsn_add_shapefile_polylines(wks,plot,shpname,plres)

;-- draw the plot and advance the frame


Document Actions