DKRZ NCL JPEG images as background map overlayed by clouds example#
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/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
Result:
