DKRZ PyNGL example overlay two different grid resolutionsΒΆ

Example:

#
#  File:
#    attach_map_to_map_Satellite.py
#
#  Synopsis:
#    Overlay different map and contour plots on a base map.
#
#  Category:
#    map plot
#    contour plot
#    polylines
#    text
#    overlay
#
#  Based on DKRZ NCL example:
#    attach_map_to_map_Satellite.ncl
#
#  Author:
#    Karin Meier-Fleischer
#
#  Date of initial publication:
#    December, 2018
#
#  Description:
#    The example shows the complexity of overlaying different maps and contour plots on a base map.
#    To point the Satellite view up a grey shadow is attached to the regional contour plot.
#
#  Effects illustrated:
#    o  Creating a map plot
#    o  Overlay two maps
#    o  Add polylines
#    o  Add text
#
#  Output:
#     A single visualization is produced.
#
'''
  PyNGL Example:     attach_map_to_map_Satellite.py

  -  Creating a map plot
  -  Overlay two maps
  -  Add polylines
  -  Add text

'''
from __future__ import print_function
import numpy as np
import os, sys
import Ngl,Nio

#--  global data
fname1 = "orog_fx_MPI-ESM-LR_rcp26_r0i0p0.nc"    #-- orography
mname1 = "sftlf_fx_MPI-ESM-LR_rcp26_r0i0p0.nc"   #-- land area fraction

#-----------------------------------------------------------------------
#-- Function: main
#-----------------------------------------------------------------------
def main():

    #--  open file and read variables
    f1    =  Nio.open_file(fname1,"r")           #-- open data file
    var1  =  f1.variables["orog"][:,:]           #-- read variable orog
    print(f1)

    mask1 =  Nio.open_file(mname1,"r")           #-- open data file
    lsm1  =  mask1.variables["sftlf"][:,:]       #-- read variable sftlf
    lsm1  =  np.where(lsm1 > 0.5, lsm1, 0)

    land_only1 = np.where(lsm1 > 0.5, var1, -101) #-- min contour value - 1

    lat   =  f1.variables["lat"][:]              #-- latitudes
    lon   =  f1.variables["lon"][:]              #-- longitudes
    dlat  =  lat[1]-lat[0]
    dlon  =  lon[1]-lon[0]

    #---------------------------------------------------------------------------
    #-- regional data
    #---------------------------------------------------------------------------
    fname2 = "/Users/k204045/data/CORDEX/EUR-11/orog_EUR-11_MPI-M-MPI-ESM-LR_historical_r0i0p0_CLMcom-CCLM4-8-17_v1_fx.nc"
    mname2 = "/Users/k204045/data/CORDEX/EUR-11/sftlf_EUR-11_MPI-M-MPI-ESM-LR_historical_r0i0p0_CLMcom-CCLM4-8-17_v1_fx.nc"

    f2     = Nio.open_file(fname2,"r")           #-- open data file
    var2   = f2.variables["orog"][:,:]           #-- read variable orog

    mask2  = Nio.open_file(mname2,"r")           #-- open data file
    lsm2   = mask2.variables["sftlf"][:,:]       #-- read variable sftlf
    lsm2   = np.where(lsm2 > 0.5, lsm2, 0)

    land_only2 = np.where(lsm2 > 0.5, var2, -101) #-- min contour value - 1

    lat2d  = f2.variables["lat"][:,:]            #-- retrieve latitudes
    lon2d  = f2.variables["lon"][:,:]            #-- retrieve longitudes

    nlat   = len(lat2d[:,0])                     #-- number of latitudes
    nlon   = len(lon2d[0,:])                     #-- number of longitudes

    #---------------------------------------------------------------------------
    #-- open workstation
    #---------------------------------------------------------------------------
    wkres                   =  Ngl.Resources()
    wkres.wkBackgroundColor = "grey20"       #-- background color
    wkres.wkForegroundColor = "white"        #-- foreground color
    wks_type = "png"                         #-- output graphic type
    wks =  Ngl.open_wks(wks_type,"plot_attach_map_to_map_Satellite",wkres) #-- open workstation

    #---------------------------------------------------------------------------
    #-- set resources
    #---------------------------------------------------------------------------
    res                             =  Ngl.Resources()
    res.nglDraw                     =  False          #-- don't draw plot yet
    res.nglFrame                    =  False          #-- don't advance frame
    res.nglMaximize                 =  False

    res.vpWidthF                    =  0.7            #-- viewport width
    res.vpHeightF                   =  0.7            #-- viewport height
    res.vpXF                        =  0.18            #-- start x-position
    res.vpYF                        =  0.85           #-- start y-position

    #---------------------------------------------------------------------------
    #-- common map resources
    #---------------------------------------------------------------------------
    mapres                          =  res

    mapres.mpDataBaseVersion        = "MediumRes"     #-- map resolution
    mapres.mpProjection             = "Satellite"     #-- map projection

    mapres.mpLimitMode              = "angles"        #-- angles represent angular
    mapres.mpLeftAngleF             =  25.            #--   deviation from the line
    mapres.mpRightAngleF            =  25.            #--   of sight from the satellite
    mapres.mpTopAngleF              =  25.            #--   to the projection center
    mapres.mpBottomAngleF           =  25.            #

    mapres.mpCenterLatF             =  40.            #-- center latitude
    mapres.mpCenterLonF             =  10.            #-- center longitude

    mapres.mpSatelliteAngle1F       =  7.*57.2957795130823*np.arcsin(1./30.)/8.    #-- 1.671437 to 3.58404
    mapres.mpSatelliteAngle2F       =  85.            #-- orthogonal view
    mapres.mpSatelliteDistF         =  2.32           #-- n-times Earth distance

    mapres.mpLabelsOn               =  False          #-- no map labels
    mapres.mpPerimOn                =  False          #-- don't draw a box around the map plot
    mapres.mpOutlineOn              =  True           #-- draw coastal outlines
    mapres.mpGeophysicalLineThicknessF = 1.7          #-- line thickness
    mapres.mpGeophysicalLineColor   = "grey30"        #-- line color
    mapres.mpGridAndLimbOn          =  False          #-- don't draw map grid lines

    mapres.mpFillOn                  =  True           #-- turn map fill on
    mapres.mpFillDrawOrder           = "PreDraw"       #-- draw filled map first
    mapres.mpOceanFillColor          = [ 0.53, 0.808, 1.0 ] #-- ocean color
    mapres.mpInlandWaterFillColor    = [ 0.53, 0.808, 1.0 ] #-- inland water color
    mapres.mpLandFillColor           = [ 0.7,  0.7,   0.7 ] #-- land color

    mapres.tmXBOn                   = False
    mapres.tmXTOn                   = False
    mapres.tmYROn                   = False
    mapres.tmYLOn                   = False
    mapres.tmXBBorderOn             = False
    mapres.tmXTBorderOn             = False
    mapres.tmYRBorderOn             = False
    mapres.tmYLBorderOn             = False

    map = Ngl.map(wks,mapres)                      #-- create the global base map
    print("---> global map done")

    #---------------------------------------------------------------------------
    #-- common contour resources
    #---------------------------------------------------------------------------
    cmap = Ngl.read_colormap_file("OceanLakeLandSnow")
    cmap[0,:] = [ 0.53, 0.808, 1.0,1.0 ]                #-- use a brighter blue

    cnres                           =  res
    cnres.cnFillOn                  =  True           #-- turn contour fill on
    cnres.cnFillMode                = "CellFill"    #-- use raster fill
    cnres.cnCellFillEdgeColor       = "gray50"

    cnres.cnLinesOn                 =  False          #-- don't draw contour lines
    cnres.cnLineLabelsOn            =  False          #-- don't draw contour line labels
    cnres.cnInfoLabelOn             =  False          #-- don't draw contour info label
    cnres.cnFillPalette             =  cmap           #-- choose color map
    cnres.cnRasterSmoothingOn       =  False          #-- do not smooth contouring
    cnres.cnLevelSelectionMode      = "ManualLevels"  #-- use manual levels
    cnres.cnMinLevelValF            =  -100.0         #-- minimum contour value
    cnres.cnMaxLevelValF            =  3000.          #-- maximum contour value
    cnres.cnLevelSpacingF           =  50.            #-- contour interval

    cnres.tiXAxisString             = ""              #-- don't draw x-axis title
    cnres.tiYAxisString             = ""              #-- don't draw y-axis title
    cnres.tiMainString              = "Overlay: regional on global grid (orography)"
    cnres.tiMainFontHeightF         =  0.02

    cnres.lbAutoManage              =  False          #-- control labelbar manually
    cnres.lbBoxMinorExtentF         =  0.15            #-- smaller labelbar boxes
    cnres.lbOrientation             = "vertical"      #-- labelbar orientation
    cnres.lbLabelFontHeightF        =  0.012          #-- labelbar font size
    cnres.pmLabelBarOrthogonalPosF  = -1.3            #-- move labelbar to the left side of plot
    cnres.pmLabelBarParallelPosF    =  0.49           #-- move labelbar downwards# default: 0.5
    cnres.lbLabelPosition           = "left"          #-- move labels to left side of labelbar
    cnres.lbRightMarginF            = -0.28           #-- move labelbar to the left

    land_only1cyc,loncyc = Ngl.add_cyclic(land_only1,lon)  #-- add cyclic point

    cnres.sfXArray                 =  loncyc
    cnres.sfYArray                 =  lat

    plot1 = Ngl.contour(wks,land_only1cyc,cnres)    #-- create global plot
    print("---> contour plot1 done")

    #-- overlay plot1 on map
    Ngl.overlay(map,plot1)
    print("---> overlay contour plot1 on global map")

    #---------------------------------------------------------------------------
    #-- regional contour resources
    #---------------------------------------------------------------------------
    cnres2                          =  cnres
    cnres2.trYReverse               =  True           #-- reverse y-axis
    cnres2.tfDoNDCOverlay           =  False          #-- transform to standard lat/lon
    cnres2.sfXArray                 =  lon2d          #-- longitude array
    cnres2.sfYArray                 =  lat2d          #-- latitude array
    cnres2.lbLabelBarOn             =  False          #-- don't draw a labelbar
    cnres2.cnFillDrawOrder          = "Draw"      #-- draw last# default: "Draw"
    cnres2.cnCellFillEdgeColor      = "transparent"

    cnres2.mpOutlineOn              =  True           #-- draw coastal outlines
    cnres2.mpGeophysicalLineThicknessF = 1.7          #-- line thickness
    cnres2.mpGeophysicalLineColor   = "gray30"        #-- line color
    cnres2.mpGridAndLimbOn          =  False          #-- don't draw map grid lines
    cnres2.mpFillOn                 =  False          #-- don't fill map

    plot2 = Ngl.contour(wks,land_only2,cnres2)    #-- create Europe plot
    print("---> regional contour plot2 done")

    #-- overlay plot2 on map
    Ngl.overlay(map,plot2)
    print("---> overlay contour plot2 on global map")

    #---------------------------------------------------------------------------
    #-- draw the edges of the regional map in black
    #---------------------------------------------------------------------------

    plres                   =  Ngl.Resources()
    plres.tfPolyDrawOrder   = "PostDraw"        #-- draw line last
    plres.gsLineThicknessF  =  1.5              #-- line thickness
    plres.gsLineColor       = "grey15"          #-- line color
    #---------------------------------------------------------------------------
    #-- edges for the shadow box and the borderline of plot2
    #---------------------------------------------------------------------------
    lon_val_lower = lon2d[0,:]
    lon_val_right = lon2d[:,nlon-1]
    lon_val_left  = lon2d[:,0]
    lon_val_upper = lon2d[nlat-1,:]

    lat_val_lower = lat2d[0,:]
    lat_val_right = lat2d[:,nlon-1]
    lat_val_left  = lat2d[:,0]
    lat_val_upper = lat2d[nlat-1,:]

    upper = Ngl.add_polyline(wks, map, lon_val_upper, lat_val_upper, plres) #-- attach upper line
    lower = Ngl.add_polyline(wks, map, lon_val_lower, lat_val_lower, plres) #-- attach lower line
    left  = Ngl.add_polyline(wks, map, lon_val_left,  lat_val_left,  plres) #-- attach left line
    right = Ngl.add_polyline(wks, map, lon_val_right, lat_val_right, plres) #-- attach right line
    print("---> black rectangle/edge done")

    #-- add units and copyrights to wks
    vpx = Ngl.get_float(map,"vpXF")             #-- retrieve value of res.vpXF from plot
    vpy = Ngl.get_float(map,"vpYF")             #-- retrieve value of res.vpYF from plot
    vpw = Ngl.get_float(map,"vpWidthF")         #-- retrieve value of res.vpWidthF from plot
    vph = Ngl.get_float(map,"vpHeightF")        #-- retrieve value of res.vpHeightF from plot

    units = f1.variables["orog"].attributes['units']

    txres = Ngl.Resources()
    txres.txFontHeightF = 0.014                 #-- font size for left, center and right string
    txres.txJust        = "CenterCenter"        #-- text justification
    Ngl.text_ndc(wks, "["+units+"]", vpx-0.04, vpy, txres)  #-- add text to wks

    txres.txFontHeightF = 0.012                 #-- font size for left, center and right string
    Ngl.text_ndc(wks,"~F35~c ~F21~~N~DKRZ", vpx+vpw-0.07, vpy-vph, txres) #-- plot copyright info

    #---------------------------------------------------------------------------
    #-- draw the map
    #---------------------------------------------------------------------------
    Ngl.draw(map)
    Ngl.frame(wks)

    Ngl.end()

if __name__ == '__main__':
    main()

Result:

PyNgl attach map to map w400