Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL JPEG images as background map overlayed by clouds example
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

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

Artikelaktionen