Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL regional map overlayed on global map Satellite projection example
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

DKRZ NCL regional map overlayed on global map Satellite projection example

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.

Example script:

;-----------------------------------------------------------
; DKRZ NCL example:  attach_map_to_map.ncl
;
; Description:      Plot a regional map with 0.11 degrees
;                   resolution on a global map with 1.875s
;                   degree resolution
;
;                    - Satellite projection
;                   - define user procedure to draw grid
;                     lines from data grid on global map
;                   - turn off grid lines on regional map
;                    - geophysical map outlines on both maps
;                   - common color map
;                    - common labelbar
;                   - draw a shadow of the regional plot
;                   - use gsn_text_ndc to draw title and
;                     copyright info strings
;
; 29.04.16  meier-fleischer(at)dkrz.de
;-----------------------------------------------------------

;-----------------------------------------------------------
;-- Procedure:  draw_data_grid_lines
;--
;--             draw grid lines of origin latitudes and
;--             longitudes
;-----------------------------------------------------------
undef("draw_data_grid_lines")
procedure draw_data_grid_lines(wks,map,var)
local lats, lons, numlon, numlat, lon_step, lat_step, ilat, ilon, x, y

begin
  BOXES = 1                                 ;-- 0: lat/lon grid lines in cell center
                                            ;-- 1: lat/lon grid lines around cells (=boxes)
  lats     =  var&lat
  lons     =  var&lon
  lat_step =  lats(1) - lats(0)
  lon_step =  lons(1) - lons(0)
  numlat   =  dimsizes(lats)
  numlon   =  dimsizes(lons)

  lonadd =  new(numlon+1,typeof(lons))      ;-- max. longitude less 360. so we have to add it
  lonadd(0:numlon-1) = lons
  lonadd(numlon)     = lons(numlon-1)+lon_step

;-- grid line resources
  gres                  =  True
  gres@gsLineColor      = "grey50"          ;-- grid line color
  gres@gsLineThicknessF =  1.7              ;-- grid line thickness

;-- assign x and x arrays for the grid lines
  x = new(numlon+1,typeof(lons))
  y = new(numlon+1,typeof(lats))

;-- draw cell centers
  if(BOXES .eq. 0) then
    dx = 0.                                 ;-- cell center
    dy = 0.                                 ;-- cell center
  else
    dx = lon_step/2.                        ;-- boxes
    dy = lat_step/2.                        ;-- boxes
  end if
 
;-- draw horizontal lines
  x = lonadd + dx
  do ilat = 0,numlat-1
       y = lats(ilat) + dy
       str = unique_string("gridline")                      ;-- create a unique name
       map@$str$ = gsn_add_polyline(wks, map, x, y, gres)   ;-- add polyline to map
  end do

;-- draw vertical lines
  y := fspan(min(lats),max(lats),numlon) + dy
  do ilon = 0,numlon
       x = lonadd(ilon) + dx
       str = unique_string("gridline")                      ;-- create unique name
       map@$str$ = gsn_add_polyline(wks, map, x, y, gres)   ;-- add polyline to map
  end do
 
end  ;-- end procedure draw_data_grid_lines


;------------
;--  Main
;------------
begin
;-- global data: open files and read data
  f1         =  addfile("$HOME/data/CMIP5/atmos/orog_fx_MPI-ESM-LR_rcp26_r0i0p0.nc","r")
  var1       =  f1->orog(:,:)
 
  mask1      =  addfile("$HOME/data/CMIP5/atmos/sftlf_fx_MPI-ESM-LR_rcp26_r0i0p0.nc","r")
  lsm1       =  mask1->sftlf(:,:)
  lsm1       =  where(lsm1.gt.0.5,-9999,lsm1)
 
  lat        =  f1->lat                             ;-- retrieve latitudes
  lon        =  f1->lon                             ;-- retrieve longitudes
  dlat       =  lat(1)-lat(0)                       ;-- lat grid interval
  dlon       =  lon(1)-lon(0)                       ;-- lon grid interval
  nlat1      =  dimsizes(lon)                       ;-- number of latitudes
  nlon1      =  dimsizes(lon)                       ;-- number of longitudes

  land_only1     =  mask(var1,lsm1,-9999)  
  land_only1!0   = "lat"                            ;-- named coordinate lat
  land_only1!1   = "lon"                            ;-- named coordinate lon
  land_only1&lat =  lat                             ;-- assign variable lat
  land_only1&lon =  lon                             ;-- assign variable lon
 
;-- regional data: open files and read data
  f2         = addfile("$HOME/data/CORDEX/EUR-11/orog_EUR-11_MPI-M-MPI-ESM-LR_historical_r0i0p0_CLMcom-CCLM4-8-17_v1_fx.nc","r")
  var2       = f2->orog(:,:)
 
  mask2      = addfile("$HOME/data/CORDEX/EUR-11/sftlf_EUR-11_MPI-M-MPI-ESM-LR_historical_r0i0p0_CLMcom-CCLM4-8-17_v1_fx.nc","r")
  lsm2       = mask2->sftlf(:,:)
  lsm2       = where(lsm2.gt.0.5,-9999,lsm2)
 
  land_only2 = var2
  land_only2 = mask(var2,lsm2,-9999)
  lat2d      = f2->lat                              ;-- retrieve latitudes
  lon2d      = f2->lon                              ;-- retrieve longitudes
  nlat       = dimsizes(lat2d(:,0))                 ;-- number of latitudes
  nlon       = dimsizes(lon2d(0,:))                 ;-- number of longitudes
 
;-- open workstation
  wks_type                        = "png"           ;-- output graphic type
  wks_type@wkBackgroundColor      = "grey20"        ;-- background color
  wks_type@wkForegroundColor      = "white"         ;-- foreground color
  wks_type@wkWidth                =  1500           ;-- workstation width
  wks_type@wkHeight               =  1500           ;-- workstation height
  wks  = gsn_open_wks(wks_type,"plot_attach_map_to_map_Satellite")

;-- common resources
  res                             =  True
  res@gsnDraw                     =  False          ;-- don't draw plot yet
  res@gsnFrame                    =  False          ;-- don't advance frame
 
;-- common map resources   
  mapres                          =  res
  mapres@vpWidthF                 =  1.0            ;-- viewport width
  mapres@vpHeightF                =  0.85           ;-- viewport height
  mapres@vpXF                     =  0.08           ;-- start x-position
  mapres@vpYF                     =  0.90           ;-- start y-position
 
  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*asin(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

;-- map resources for global map
  mpres = mapres
  mpres@mpFillOn                  =  True           ;-- turn map fill on
  mpres@mpFillDrawOrder           = "PreDraw"       ;-- draw filled map first
  mpres@mpOceanFillColor          = (/ 0.824, 0.961, 1.0 /) ;-- ocean color
  mpres@mpInlandWaterFillColor    = (/ 0.824, 0.961, 1.0 /) ;-- inland water color
  mpres@mpLandFillColor           = (/ 0.7, 0.7, 0.7 /)     ;-- land color

  map = gsn_csm_map(wks,mpres)                      ;-- create the global base map
  print("---> global map done")

;-- common contour resources
  cnres = res
  cnres@cnFillOn                  =  True           ;-- turn contour fill on
  cnres@cnFillMode                = "RasterFill"    ;-- use raster fill
  cnres@cnMissingValFillColor     = "steelblue1"    ;-- missing value color
  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             = "OceanLakeLandSnow" ;-- 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@gsnLeftString             = "["+var1@units+"]"  ;-- draw units above the labelbar
  cnres@gsnLeftStringOrthogonalPosF = -0.01         ;-- move the left string downward
  cnres@gsnLeftStringParallelPosF = -0.066          ;-- move left string left
  cnres@gsnLeftStringFontHeightF  =  0.012          ;-- left string font size
  cnres@gsnRightString            = ""              ;-- don't draw right string
  cnres@tiXAxisString             = ""              ;-- don't draw x-axis title
  cnres@tiYAxisString             = ""              ;-- don't draw y-axis title
 
  cnres@lbAutoManage              =  False          ;-- control labelbar manually
  cnres@lbBoxMinorExtentF         =  0.2            ;-- smaller labelbar boxes
  cnres@lbOrientation             = "vertical"      ;-- labelbar orientation
  cnres@lbLabelFontHeightF        =  0.012          ;-- labelbar font size
  cnres@pmLabelBarOrthogonalPosF  = -1.225          ;-- 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
 
;-- contour resources plot1
  cnres1 = cnres
  cnres1@gsnAddCyclic             =  True           ;-- add cyclic point
 
  plot1 = gsn_csm_contour(wks,land_only1,cnres1)    ;-- create global plot
  print("---> contour plot1 done")

;-- overlay plot1 on map
  overlay(map,plot1)
  print("---> overlay contour plot1 on global map")
 
;-- add grid lines from original grid
  draw_data_grid_lines(wks,map,land_only1)          ;-- procedure defined above
  print("---> draw grid lines 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          = "PostDraw"      ;-- draw last; default: "Draw"
 
  plot2 = gsn_csm_contour(wks,land_only2,cnres2)    ;-- create Europe plot
  print("---> regional contour plot2 done")

;-- edges for the shadow box and the borderline of plot2
  lon_val_upper  =  lon2d(nlat-1,:)
  lat_val_upper  =  lat2d(nlat-1,:)
  lon_val_lower  =  lon2d(0,:)
  lat_val_lower  =  lat2d(0,:)
  lon_val_left   =  lon2d(:,0)
  lat_val_left   =  lat2d(:,0)
  lon_val_right  =  lon2d(:,nlon-1)
  lat_val_right  =  lat2d(:,nlon-1)

;-- draw a "shadow" rectangle around the regional map for better visibility
  dx = 1.5                                          ;-- x-offset of shadow rectangle
  dy = 1.0                                          ;-- y-offset of shadow rectangle
 
;-- x-array for shadow rectangle
  x  = array_append_record((/lon_val_lower+dx/), (/lon_val_right+dx/),       0)
  x := array_append_record((/x/),                (/lon_val_upper(::-1)+dx/), 0)
  x := array_append_record((/x/),                (/lon_val_left(::-1)+dx/),  0)
  x := array_append_record((/x/),                (/lon2d(0,0)+dx/),          0)

;-- y-array for shadow rectangle
  y  = array_append_record((/lat_val_lower-dy/), (/lat_val_right-dy/),       0)
  y := array_append_record((/y/),                (/lat_val_upper(::-1)-dy/), 0)
  y := array_append_record((/y/),                (/lat_val_left(::-1)-dy/),  0)
  y := array_append_record((/y/),                (/lat2d(0,0)-dy/),          0)

;-- draw the transparent shadow rectangle (filled polygon)
  pgres                 =  True
  pgres@tfPolyDrawOrder = "PostDraw"                ;-- draw the rectangle last
  pgres@gsLineColor     = "grey30"                  ;-- rectangle line color
  pgres@gsFillColor     = "grey30"                  ;-- rectangle fill color
  pgres@gsFillOpacityF  =  0.7                      ;-- light transparency

  shadow = gsn_add_polygon(wks, map, x, y, pgres)   ;-- attach shadow rectangle to map
  print("---> regional plot2 shadow done")
 
;-- overlay regional plot on map
  overlay(map,plot2)
  print("---> overlay regional plot2 done")

;-- create a simple global map with geophysical outlines only
  mpres2 = mapres
  mpres2@mpGridAndLimbOn  =  False                  ;-- don't draw map grid lines
  mpres2@mpFillOn         =  False                  ;-- don't fill map
 
  map2 = gsn_csm_map(wks,mpres2)                    ;-- create simple map
  print("---> simple map2 outlines done")
 
;-- attach map2 to map to draw the geophysical outlines ontop of the regional map, too.                                   
  anres =  True
  anid  =  gsn_add_annotation(map,map2,anres)       ;-- attach map2 to map
  print("---> add map2 annotation done")

;-- draw the edges of the regional map in black
  plres                   =  True
  plres@tfPolyDrawOrder   = "PostDraw"              ;-- draw line last
  plres@gsLineThicknessF  =  1.5                    ;-- line thickness
  plres@gsLineColor       = "grey15"                ;-- line color

  upper = gsn_add_polyline(wks, map2, lon_val_upper, lat_val_upper, plres) ;-- attach upper line
  lower = gsn_add_polyline(wks, map2, lon_val_lower, lat_val_lower, plres) ;-- attach lower line
  left  = gsn_add_polyline(wks, map2, lon_val_left,  lat_val_left,  plres) ;-- attach left line
  right = gsn_add_polyline(wks, map2, lon_val_right, lat_val_right, plres) ;-- attach right line
  print("---> black rectangle/edge done")

;-- draw the map
  draw(map)

;-- draw manual title on top of the plot
  title = "Overlay: regional on global grid (orography)" ;-- title string
  tires                   =  True                   ;-- text resources copyright string
  tires@txFontColor       = "white"                 ;-- change to white
  tires@txJust            = "CenterCenter"          ;-- text justification
  tires@txFontHeightF     =  0.024                  ;-- make font size smaller
  tires@txFontThicknessF  =  3.0                    ;-- thicker font
  gsn_text_ndc(wks,title, 0.6, 0.95, tires)         ;-- plot copyright info

;-- add copyright info string
  txres                   =  True                   ;-- text resources copyright string
  txres@txFontColor       = "white"                 ;-- change to white
  txres@txJust            = "BottomRight"           ;-- text justification
  txres@txFontHeightF     =  0.013                  ;-- make font size smaller
  txres@txFontThicknessF  =  3.0                    ;-- thicker font
  gsn_text_ndc(wks,"~F35~c ~F21~~N~DKRZ", 0.98, 0.013, txres) ;-- plot copyright info
 
;-- advance the frame
  frame(wks)
  print("---> all done")

end

Artikelaktionen