You are here: Home / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL JPEG images as background map overlayed by clouds example

DKRZ NCL JPEG images as background map overlayed by clouds example

This example demonstrates how to read in a JPEG file which is used as background map. The ICON cloud data are drawn on top of the map using different color transparencies.

Example script:

;-------------------------------------------------------------------
;-- DKRZ NCL example:
;--   NCL_ICON_triangles_cloud_cover_earth_black_white_trans_opaq.ncl  
;--
;--   Read ICON data and plot variable 'clt'.
;--
;--   Input data file does not contain the grid information, a second
;--   file r2b4_amip.nc with the grid information must be opened.
;--   Original grid files on blizzard:
;--              /pool/data/ICON/grids/grids/private/r2b4_amip
;--
;--   Script based on the NCL example icon_5.ncl from
;--           http://ncl.ucar.edu/Applications/Scripts/icon_5.ncl
;--
;--   modified:
;--      - don't wrap lines around for 'low' and 'high' levels as well
;--           (original script regards only to 'middle' levels which  
;--            caused polygon line wrapping)
;--      - optimize the loop for 'low', 'high' and 'middle' levels
;--      - considering missing values (here: filled with 'red' color)
;--      - underlaying JPG image of the Earth
;--      - turn off/on drawing the edges
;--      - use transparent to opaque white colormap for variable 'clt'
;--
;-- 03.01.14  meier-fleischer(at)dkrz.de
;-------------------------------------------------------------------
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

;----------------------------------------------------------------------
;-- set some 'global' variables
;----------------------------------------------------------------------
plotfile = "plot_ICON_triangles_cloud_cover_earth_black_white_trans_opaq"


;----------------------------------------------------------------------
;-- Function:  recreate_jpeg_image
;--
;-- Recreate_jpeg_image imports a JPEG image  of the Earth and recreate
;-- an NCL object file for plotting.
;-- 
;-- Based on NCL example topo_9.ncl
;----------------------------------------------------------------------
undef("recreate_jpeg_image")
function recreate_jpeg_image(wks,var,minlat,maxlat,minlon,maxlon)
local red, greens, blues, ramp, band_dims, nlon, nlat, Band1, Band2, Band3
begin
  print("-- procedure recreate_jpeg_image")

  nc_filename = "$HOME/data/EarthMap_2500x1250.nc"
  imgfile     = addfile(nc_filename,"r")

;-- Read the three bands of data
  Band1     = where(imgfile->Band1.gt.255, 255, imgfile->Band1)  ; red channel
  Band2     = where(imgfile->Band2.gt.255, 255, imgfile->Band2)  ; green channel
  Band3     = where(imgfile->Band3.gt.255, 255, imgfile->Band3)  ; blue channel
  band_dims = dimsizes(Band3)
  nlat      = band_dims(0)
  nlon      = band_dims(1)
  print("     dimensions of Earth image: " + nlat + " x " + nlon)

;-- add lat/lon data so we can overlay on a map, and/or overlay contours. 
;-- We know the image is global, cylindrical equidistant, and centered 
;-- about lon=0.
  lat       =  fspan( -90, 90,nlat)                    ;-- generate lat array
  lon       =  fspan(-180,180,nlon)                    ;-- generate lon array
  lon!0     = "lon"                                    ;-- named dimension lon
  lat!0     = "lat"                                    ;-- named dimension lat
  lon@units = "degrees_east"                           ;-- lon units
  lat@units = "degrees_north"                          ;-- lat units

  Band1!0   = "lat"                                    ;-- named dimension lat
  Band1!1   = "lon"                                    ;-- named dimension lon
  Band2!0   = "lat"
  Band2!1   = "lon"
  Band3!0   = "lat"
  Band3!1   = "lon"
  Band1&lat =  lat                                     ;-- set lat data
  Band1&lon =  lon                                     ;-- set lon data
  Band2&lat =  lat
  Band2&lon =  lon
  Band3&lat =  lat
  Band3&lon =  lon
  
;-- set plotting resources
  eres                       =  True
  eres@gsnDraw               =  False                  ;-- don't draw plot yet
  eres@gsnFrame              =  False                  ;-- don't advance the frame
  eres@gsnMaximize           =  True                   ;-- maximize the plot output
  eres@gsnLeftString         = ""                      ;-- don't draw left string
  eres@gsnRightString        = ""                      ;-- don't draw right string

  eres@cnFillOn              =  True                   ;-- use filled contour
  eres@cnFillMode            = "RasterFill"            ;-- use RasterFill, its faster
  eres@cnLevelSelectionMode  = "EqualSpacedLevels"     ;-- equal spaced contour levels
  eres@cnMaxLevelCount       =  254                    ;-- max number of contour levels
  eres@cnFillBackgroundColor =  (/ 1., 1., 1., 1./)    ;-- use background color
  eres@cnLinesOn             =  False                  ;-- draw contour lines
  eres@cnLineLabelsOn        =  False                  ;-- don't draw contour labels
  eres@cnInfoLabelOn         =  False                  ;-- don't draw info label

  eres@lbLabelBarOn          =  False                  ;-- don't draw a labelbar for the image
 
;-- construct RGBA colormaps...
  ramp        = fspan(0., 1., 255)                     ;-- define 255 incremental steps
  reds        = new((/255, 4/), float)                 ;-- assign reds color array
  greens      = new((/255, 4/), float)                 ;-- assign greens color array
  blues       = new((/255, 4/), float)                 ;-- assign blues color array
  reds        = 0                                      ;-- initialize all reds elements to 0
  greens      = 0                                      ;-- initialize all greens elements to 0
  blues       = 0                                      ;-- initialize all blues elements to 0
  reds(:,0)   = ramp                                   ;-- RGB, red=ramp, green=0,    blue=0
  greens(:,1) = ramp                                   ;-- RGB, red=0,    green=ramp, blue=0
  blues(:,2)  = ramp                                   ;-- RGB, red=0,    green=0,    blue=ramp

;-- the red contour map is plotted fully opaque; the green and blue are plotted 
;-- completely transparent.
;-- When overlain, the colors combine (rather magically).
  reds(:,3)   = 1.                                     ;-- RGBA  A=1. (100% opaque)
  greens(:,3) = 0                                      ;-- RGBA  A=0  (100% transparent)
  blues(:,3)  = 0                                      ;-- RGBA  A=0  (100% transparent)

  eres@cnFillColors = greens 
  greenMap = gsn_csm_contour(wks, Band2, eres)         ;-- plot all greens values, no map

  eres@cnFillColors = blues
  blueMap = gsn_csm_contour(wks, Band3, eres)          ;-- plot all blues values, no map
 
;-- this will be our base, so make it a map plot.
  eres@cnFillColors             =  reds                ;-- use given color array
  eres@gsnAddCyclic             =  False               ;-- don't add cyclic point
  eres@gsnLeftString            =  var@long_name       ;-- use var long_name
  eres@gsnRightString           =  var@units           ;-- use var units
  eres@mpFillOn                 =  False               ;-- don't draw filled land
  eres@mpDataBaseVersion        = "MediumRes"          ;-- map data base
  eres@mpLimitMode              = "LatLon"             ;-- map limits mode
  eres@mpMinLatF                =  minlat              ;-- lon minimum
  eres@mpMaxLatF                =  maxlat              ;-- lon maximum
  eres@mpMinLonF                =  minlon              ;-- lat minimum
  eres@mpMaxLonF                =  maxlon              ;-- lat maximum
  eres@mpCenterLonF             =  0                   ;-- map center lon
  eres@mpGridLonSpacingF        =  30.                 ;-- x-tickmark grid spacing
  eres@mpGridLatSpacingF        =  30.                 ;-- y-tickmark grid spacing
  eres@tmXBLabelFontHeightF     =  0.01                ;-- decrease x label font size
  eres@tmYLLabelFontHeightF     =  0.01                ;-- decrease y label font size

  redMap = gsn_csm_contour_map(wks, Band1, eres)       ;-- plot all reds values, with map

;-- overlay everything to create the topo map
  overlay(redMap, greenMap)                            ;-- overlay greenMap on redMap
  overlay(redMap, blueMap)                             ;-- overlay blueMap on redMap

  return(redMap)
end

;-----------------------------------------------------------------------
;-- Procedure:  add_lab_transparent
;-- 
;-- Create and add labelbar; solid fill with pattern underneath
;-----------------------------------------------------------------------
undef("add_lab_transparent")
procedure add_lab_transparent(wks,map,cmap,levels)
begin
  print("-- procedure add_lab_transparent")
   
  getvalues map@contour
     "vpWidthF"  : vpw                                 ;-- retrieve viewport width
     "vpHeightF" : vph                                 ;-- retrieve viewport height
     "vpXF"      : vpx                                 ;-- retrieve viewport start x
     "vpYF"      : vpy                                 ;-- retrieve viewport start y
  end getvalues
         
;-- set the labels
  labs     =  levels+""
  nlevs    =  dimsizes(labs)

;-- set labelbar width, height and position
  lbwidth  =  0.7
  lbheight =  0.08
  lbx      = ((1.0-(lbwidth+vpx))/2) + vpx/2
  lby      = ((1.0-(vpy-vph))/2) - 0.13

;-- set global labelbar resources
  lbres                       =  True
  lbres@gsnFrame              =  False                 ;-- don't advance frame
  lbres@vpWidthF              =  lbwidth               ;-- width of labelbar
  lbres@vpHeightF             =  lbheight              ;-- height of labelbar
  lbres@vpXF                  =  lbx                   ;-- labelbar x-position
  lbres@vpYF                  =  lby                   ;-- labelbar y-position
  lbres@lbPerimOn             =  False                 ;-- no label bar box
  lbres@lbOrientation         =  "Horizontal"          ;-- orientation
  lbres@lbLabelFontHeightF    =  0.015                 ;-- label font height
  lbres@lbLabelAlignment      =  "InteriorEdges"       ;-- where to label
  lbres@lbMonoFillPattern     =  True                  ;-- fill solid

;-- set pattern labelbar resources and create the labelbar and annotation
  lbres2                      =  lbres
  lbres2@lbFillPattern        =  6                     ;-- for transparent colors lay pattern 6 underneath
  lbres2@lbFillDotSizeF       =  0.01                  ;-- increase the dot size of the pattern 
  lbres2@lbFillLineThicknessF =  2.0                   ;-- increase the line thickness
  lbres2@lbFillScaleF         =  0.7                   ;-- increase the fill scale of the pattern
  lbres2@lbFillColors         = "black"                ;-- use colors
  lbres2@lbMonoFillColor      =  True                  ;-- mono fill color black
  gsn_labelbar_ndc(wks,nlevs+1,labs,lbx,lby,lbres2)    ;-- draw labelbar with black pattern

;-- set transparency labelbar resources and create the labelbar and annotation
  lbres1                      =  lbres
  lbres1@lbFillColors         =  cmap                  ;-- use colors
  lbres1@lbFillPattern        =  0                     ;-- fill solid
  lbres1@lbMonoFillColor      =  False                 ;-- no mono fill color
  gsn_labelbar_ndc(wks,nlevs+1,labs,lbx,lby,lbres1)    ;-- draw transparent labelbar

;-- set text resources and draw text on top of the labelbar
  txres                       =  True                  ;-- text resources
  txres@txFontHeightF         =  0.010                 ;-- text font size
  xl = lbx + 0.04                                      ;-- x position
  yb = lby - 0.01                                      ;-- y position
  gsn_text_ndc(wks,"opacity 0%",xl,yb, txres)          ;-- draw left string
  xr = xl + lbwidth - 0.085                            ;-- x position
  gsn_text_ndc(wks,"opacity 100%",xr,yb, txres)        ;-- draw right string
  
end

;-----------------------------------------------------------------------
;-- Procedure:  plot_ICON_triangles
;-- 
;-- Plot ICON data as triangles or contour lines
;-----------------------------------------------------------------------
undef("plot_ICON_triangles")
procedure plot_ICON_triangles(wks,var,x,y,vlon,vlat,res_in,cmap)
local plot, cnres, mpres, pres
begin
  print("-- procedure plot_ICON_triangles")
  nv        =  dimsizes(vlon(0,:))                ;-- no of points in polygon
  lonMin    = -180                                ;-- longitude minimum
  lonMax    =  180                                ;-- longitude maximum
  latMin    =  -90                                ;-- latitude minimum
  latMax    =   90                                ;-- latitude maximum
  mapCenter =    0                                ;-- center of map

  print("     Data longitude min/max: " + min(vlon) + "   " + max(vlon))
  print("     Cell points:            " + nv)
  print("     Plot area:              "+lonMin+","+lonMax+" "+latMin+","+latMax)

;-- set contour resources
  cnres                      =  res_in
  cnres@gsnDraw              =  False             ;-- don't draw the plot
  cnres@gsnFrame             =  False             ;-- don't advance the frame
  cnres@gsnMaximize          =  True              ;-- maximize plot output
    
  cnres@sfXArray             =  x                 ;-- transform x to mesh scalar field
  cnres@sfYArray             =  y                 ;-- transform y to mesh scalar field

  cnres@cnFillPalette        =  cmap              ;-- use cmap colormap
  cnres@cnLinesOn            =  False             ;-- don't draw contour lines
  cnres@cnInfoLabelOn        =  False             ;-- switch off contour info label
  cnres@cnFillOn             =  True              ;-- contour fill on
  cnres@cnFillMode           = "RasterFill"       ;-- set fill mode
  cnres@cnRasterSmoothingOn  =  True              ;-- smooth contouring
  cnres@cnMissingValFillColor = "red"             ;-- draw missing values in red
  cnres@lbLabelBarOn         =  False

  plot = gsn_csm_contour(wks,var,cnres)           ;-- create the contour plot, but don't draw it

;-- to plot the filled polygons we need the data levels and their colors
  getvalues plot
      "cnLevels"     :  levels                    ;-- retrieve the levels from plot
      "cnFillColors" :  colors                    ;-- retrieve the colors from plot
  end getvalues
  nlevels = dimsizes(levels)                      ;-- number of levels

  plot = setColorContourClear(plot,min(var),max(var)) ;-- clear plot, but keep all the information

;-- set map resources
  mpres                      =  True
  mpres@gsnDraw              =  False             ;-- don't draw the plot
  mpres@gsnFrame             =  False             ;-- don't advance the frame
  mpres@gsnMaximize          =  True              ;-- maximize plot output
  mpres@gsnPaperOrientation  = "landscape"
  mpres@gsnLeftString        = " " 

  mpres@mpLimitMode          = "LatLon"           ;-- map limits mode
  mpres@mpMinLatF            =  latMin            ;-- latitude minimum
  mpres@mpMaxLatF            =  latMax            ;-- latitude maximum
  mpres@mpMinLonF            =  lonMin            ;-- longitude minimum
  mpres@mpMaxLonF            =  lonMax            ;-- longitude maximum
  mpres@mpCenterLonF         =  0                 ;-- map center lon
  mpres@mpFillOn             =  False             ;-- don't draw filled map
  mpres@mpGridLonSpacingF    =  30.
  mpres@mpGridLatSpacingF    =  30.

;-- create map using the JPEG file of the Earth
  map = recreate_jpeg_image(wks,var,latMin,latMax,lonMin,lonMax)

;-- add labelbar using transparency/opaque color with underlaying dashpattern
  add_lab_transparent(wks,map,colors,levels+"")

;-- create color array for triangles
  ntri_calc = 0
  ntri      = dimsizes(y)                         ;-- number of triangles
  gscolors  = new(ntri,integer)
  gscolors  = -1                                  ;-- initialize to transparent

;-- set resources for the triangles (polygons)
  pres                       =  True
  pres@gsEdgesOn             =  False             ;-- turn on edges
  pres@gsFillIndex           =  0                 ;-- solid fill

;-- set color for data less than given minimum value var_min
  vlow = ind(var .lt. levels(0))                  ;-- get the indices of values less levels(0)
  if (.not. ismissing(vlow)) then
     gscolors(vlow) = colors(0)                   ;-- choose color
     ntri_calc = dimsizes(vlow)                   ;-- number of triangles
  end if
;-- set colors for all cells in between var_min and var_max
  do i = 1, dimsizes(levels) - 1
     vind := ind(var .ge. levels(i-1) .and. var .lt. levels(i))  ;-- get the indices of 'middle' values
     gscolors(vind) = colors(i)                   ;-- choose the colors
     ntri_calc = ntri_calc + dimsizes(vind)       ;-- number of triangles
  end do

;-- set color for data greater than given maximum var_max
  nc=dimsizes(colors)-1                           ;-- get the number of colors minus one
  nl=dimsizes(levels)-1                           ;-- get the number of levels minues one
  vhgh := ind(var .gt. levels(nl))                ;-- get indices of values greater levels(nl)
  gscolors(vhgh) = colors(nc)                     ;-- choose color
  ntri_calc = ntri_calc + dimsizes(vhgh)          ;-- number of triangles

  print("--> calculated cells:  "+ ntri_calc)

;-- Attach all the triangles using the list of colors
  pres@gsColors   = gscolors
  pres@gsSegments = ispan(0,dimsizes(var) * 3,3)  ;-- assign segments array
  polygon = gsn_add_polygon(wks,map,ndtooned(vlon),ndtooned(vlat),pres)  ;-- draw all triangles

  draw(map)                                          ;-- draw the map

end

;=======================================================================
;-- main script
;=======================================================================
begin
  DataFileName = "$HOME/data/ICON/atm_phy_mag0004_1985.nc"  ;-- data file
  gridinfofile = "$HOME/data/ICON/grids/r2b4_amip.nc"       ;-- grid info file

  title        = "ICON:  overlay on satellite image of the earth" ;-- title string
  bottom_string = "$HOME/NCL_support/ICON_triangles_cloud_cover_earth_black_white_trans_opaq.ncl" ;-- bottom info
  
  VarName      = "clt"                               ;-- variable name       
  scale        =  100                                ;-- multiply by

  misscol      = "red"                               ;-- color for missing values
  
  varMin       =    0                                ;-- minimum contour level
  varMax       =   90                                ;-- maximum contour level
  varInt       =    5                                ;-- interval between contours

  lonMin       = -180                                ;-- longitude minimum
  lonMax       =  180                                ;-- longitude maximum
  latMin       =  -90                                ;-- latitude minimum
  latMax       =   90                                ;-- latitude maximum
  
;-------------------------------------------------------------------
;-- add data and grid file
;-------------------------------------------------------------------
  File         = addfile( DataFileName, "r" )        ;-- add data file
  GridInfoFile = addfile( gridinfofile, "r" )        ;-- add grid file (not contained in data file!!!)
  var          = File->$VarName$(time|0,ncells|:)    ;-- set variable with dims: (time,ncell)
  var          = var*scale                           ;-- scale factor (here 100 to get %)
  var@units    = "%"                                 ;-- set new variable units

;-------------------------------------------------------------------
;-- define the x-, y-values and the polygon points
;-------------------------------------------------------------------
  rad2deg =  45./atan(1.)                            ;-- radians to degrees
 
  x       = GridInfoFile->clon *rad2deg              ;-- cell center, lon
  y       = GridInfoFile->clat *rad2deg              ;-- cell center, lat
  x!0     = "lon"                                    ;-- set named dimension lon
  y!0     = "lat"                                    ;-- set named dimension lat
  x@units = "degrees_east"                           ;-- set lon units
  y@units = "degrees_north"                          ;-- set lat units

  vlon    = GridInfoFile->clon_vertices * rad2deg    ;-- cell longitude vertices
  vlon    = where(vlon.lt.0, vlon + 360, vlon)       ;-- longitude: 0-360
  vlat    = GridInfoFile->clat_vertices * rad2deg    ;-- cell lattitude vertices
  nv      = dimsizes(vlon(0,:))                      ;-- number of points in polygon

  print("")
  print("Data longitude min/max: " + min(vlon) + "   " + max(vlon) + "   cell points = " + nv)
  print("")
  print("Plot area: "+lonMin+","+lonMax+" "+latMin+","+latMax)
  print("")
  print("Min: "+min(var)+"   Max: "+max(var))
  print("")
 
;-- generate the levels and labels
  ncols   = toint(((varMax-varMin)/varInt)+2)        ;-- number of colors
  npts    = ((varMax-varMin)/varInt)+1               ;-- number of points
  levels  = fspan(varMin,varMax,npts)                ;-- assign levels array
  labels  = levels+""                                ;-- assign labels array
  nlevels = dimsizes(levels)                         ;-- number of levels

;-- generate white color map with transparency (RGBA)
  cmap      = new((/ncols,4/), "float")              ;-- assign colormap array
  cmap      = 1.                                     ;-- set all colors to white and 100% opaque
  do i=0,ncols-1
     cmap(i,3) = i*(1./(ncols-1))                    ;-- increase opacity
  end do

;-------------------------------------------------------------------
;-- open workstation
;-------------------------------------------------------------------
  wks_type            = "png"
  wks_type@wkWidth    =  1024                        ;-- set workstation width in pixel
  wks_type@wkHeight   =  1024                        ;-- set workstation height in pixel
  wks_type@wkBackgroundColor = "gray70"
  wks_type@wkForegroundColor = "black"
  wks = gsn_open_wks(wks_type,plotfile)

;-------------------------------------------------------------------
;-- set resources to get the levels and colors for the triangles
;-------------------------------------------------------------------
  res                      =  True
  res@cnLevelSelectionMode = "ManualLevels"          ;-- set manual contouring levels on
  res@cnLevelSpacingF      =  varInt                 ;-- contouring interval
  res@cnMinLevelValF       =  varMin                 ;-- contouring minimum
  res@cnMaxLevelValF       =  varMax                 ;-- contouring maximum

  plot_ICON_triangles(wks,var,x,y,vlon,vlat,res,cmap);-- plot the ICON data as triangles

;-------------------------------------------------------------------
;-- write title on top of the plot
;-------------------------------------------------------------------
  tires                =  True                       ;-- text resources
  tires@txFontHeightF  =  0.015                      ;-- text font size
  tires@txJust         = "CenterCenter"              ;-- text justification
  tires@txFont         = "courier"                   ;-- text font
  xl = 0.5                                           ;-- text x position
  yb = 0.78                                          ;-- text y position
  gsn_text_ndc(wks,title,xl,yb, tires)               ;-- plot title string on top

;-------------------------------------------------------------------
;-- just hold a small border around the plot (later we will cut off not used space from plot)
;-------------------------------------------------------------------
  plres                  =  True
  plres@gsLineColor      = "black"                 ;-- color of lines
  plres@gsLineThicknessF =  2.0                    ;-- thickness of lines
  
  plx = (/0.016,0.99,0.99,0.016,0.016/)
  ply = (/0.15,0.15,0.83,0.83,0.15/)
  gsn_polyline_ndc(wks,plx,ply,plres)
  
  frame(wks)                                         ;-- advance the frame
;-------------------------------------------------------------------

;-------------------------------------------------------------------
;-- cut off the  white border
;-------------------------------------------------------------------
  cmd = "convert -alpha off -density 300 -trim "+plotfile+".png tmp9999.png"
  system(cmd)
  system("mv tmp9999.png "+plotfile+".png")
  
end

Document Actions