Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL shapefile mean temperature change German coast example
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

DKRZ NCL shapefile mean temperature change German coast example

NCL script to create a plot of averaged dummy data using shapefile containing the German coastal counties.

DKRZ NCL script:

;----------------------------------------------------------------------
; DKRZ NCL Examples:   NCL_shapefile_mean_temp_change_German_coast.ncl
;   
; Description:       Average "dummy" data over all counties of a selected 
;                    state of Germany and plot them on a map of Germany.
; 
; Library file:      $HOME/NCL/lib/lib_shapefiles.ncl
; 
; Usage:             ncl <arguments> DEU_adm3_avg_over_counties_plus_borderlines.ncl
; 
;     Arguments:     state_name=<string>                     ;--  default: "Schleswig-Holstein"
;                    county_border=<True or False>           ;--  default: True
;                    states_border=<True or False>           ;--  default: True
;                    country_border=<True or False>          ;--  default: True
;                    subregion="minlon,maxlon,minlat,maxlat" ;--  default: no sub-region
; Examples:
; 
;    1. Draw "Schleswig-Holstein (default) and plot only the sub-region
;
;         ncl 'subregion="7.8,11.9,53.0,55.3"' DEU_adm3_avg_over_ALL_counties_plus_borderlines.ncl
; 
;    2. Draw "Schleswig-Holstein (default) but don't draw the border of all states
; 
;         ncl 'states_border=False' DEU_adm3_avg_over_ALL_counties_plus_borderlines.ncl
; 
;    3. Select the state "Hessen" but don't draw the borderline of Germany
; 
;         ncl 'state_name="Hessen"' 'country_border=False' DEU_adm3_avg_over_ALL_counties_plus_borderlines.ncl
;  
;    4. Coastal region
;         ncl 'subregion="6.5,14.75,50.,55.5"' DEU_adm3_avg_over_ALL_counties_plus_borderlines.ncl
; 
; 26.02.14 kmf
;----------------------------------------------------------------------
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"

;----------------------------------------------------------------------
;-- Function:  add_shapefile_polygons(...)
;--                ->  add polygons of the selected segments of the 
;--                    shapefile to the plot
;----------------------------------------------------------------------
undef("add_shapefile_polygons")
;----------------------------------------------------------------------
function add_shapefile_polygons(wks,plot,state_name,colors,f,gnres)
;----------------------------------------------------------------------
local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
      geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, i, lat, lon, \
      startSegment, numSegments, seg, startPT, endPT, npoly, npl
begin
;-- some error handling
  if(ismissing(f)) then
     print("Error: add_shapefile_polygons: Can't open shapefile '" + fname + "'")
     return(new(1,graphic))
  end if
  if([Email protection active, please enable JavaScript.]."polygon") then
     print("Error: add_shapefile_polygons: Attribute geometry_type must be 'polygon'")
     return(new(1,graphic))
  end if
  
;-- read shapefile data
  segments      =  f->segments
  geometry      =  f->geometry
  segsDims      =  dimsizes(segments)
  geomDims      =  dimsizes(geometry)
  lon           =  f->x                                   ;-- longitudes array of counties
  lat           =  f->y                                   ;-- latitudes array of counties
  geom_segIndex =  f@geom_segIndex
  geom_numSegs  =  f@geom_numSegs
  segs_xyzIndex =  f@segs_xyzIndex
  segs_numPnts  =  f@segs_numPnts
  numFeatures   =  geomDims(0)
  
;-- create array to hold all polylines
  npoly =  sum(geometry(:,geom_numSegs))                 ;-- sum of all counties polygons
  poly  =  new(npoly,graphic)                            ;-- array of all counties polygons
  npl   =  0                                             ;-- polyline counter
  j     =  0                                             ;-- counter for the colors array
  
;-- draw the color filled polygons
  do i=0, numFeatures-1  
     startSegment = geometry(i, geom_segIndex)
     numSegments  = geometry(i, geom_numSegs)
     do seg=startSegment, startSegment+numSegments-1
        startPT = segments(seg, segs_xyzIndex)
        endPT   = startPT + segments(seg, segs_numPnts) - 1
        gnres@gsFillColor     = colors(j)                 ;-- set color
        gnres@tfPolyDrawOrder = "PostDraw"
        poly(npl)  = gsn_add_polygon (wks, plot, lon(startPT:endPT),lat(startPT:endPT), gnres)
        npl = npl + 1
     end do
     j=j+1
  end do
  
  return(poly)            ;-- return the polygon array
end

;----------------------------------------------------------------------
;-- Function:  add_shapefile_polylines(...)
;--                ->  add polylines of the selected segments of the
;--                    shapefile to the plot
;----------------------------------------------------------------------
undef("add_shapefile_polylines")
;----------------------------------------------------------------------
function add_shapefile_polylines(wks,plot,state_name,colors,f,gnres)
;----------------------------------------------------------------------
local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
      geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, i, lat, lon, \
      startSegment, numSegments, seg, startPT, endPT, npoly, npl
begin
;-- some error handling
  if(ismissing(f)) then
     print("Error: add_shapefile_polys: Can't open shapefile '" + fname + "'")
     return(new(1,graphic))
  end if
  if([Email protection active, please enable JavaScript.]."polygon") then
     print("Error: add_shapefile_polys: Attribute geometry_type must be 'polygon'")
     return(new(1,graphic))
  end if
  
;-- read shapefile data
  segments      =  f->segments
  geometry      =  f->geometry
  segsDims      =  dimsizes(segments)
  geomDims      =  dimsizes(geometry)
  lon           =  f->x                                   ;-- longitudes array of counties
  lat           =  f->y                                   ;-- latitudes array of counties
  geom_segIndex =  f@geom_segIndex
  geom_numSegs  =  f@geom_numSegs
  segs_xyzIndex =  f@segs_xyzIndex
  segs_numPnts  =  f@segs_numPnts
  numFeatures   =  geomDims(0)

;-- grab the indices containing the counties of the selected state
  states        =  f->NAME_1
  names3        =  f->NAME_3
  DEU_counties  =  ind(names3.ne."")

;-- get state_name counties
  wc=new(dimsizes(names3),typeof(names3))
  if(.not.isatt(wc,"_FillValue")) then
    wc@_FillValue = default_fillvalue(typeof(names3))          ;-- make sure "wc" has a missing value
  end if
  n=0
  do m=0,dimsizes(names3)-1
    if(states(m).eq.state_name) then
       wc(n) = names3(m)                                       ;-- get all relevant counties
    else
       wc(n) = default_fillvalue(typeof(names3))               ;-- set all other to missing value
    end if
    n=n+1
  end do
  wcounties = ind(.not. ismissing(wc))
  state_counties = new(dimsizes(wcounties),string)
  do jj=0,dimsizes(wcounties)-1
     state_counties(jj) = names3(wcounties(jj))                ;-- get the names of the relevant counties
  end do

;-- create array to hold all polylines
  npoly   =  sum(geometry(:,geom_numSegs))               ;-- sum of all counties polygons
  poly    =  new(npoly,graphic)                          ;-- array of all counties polygons
  npl     =  0                                           ;-- polyline counter
  j       =  0                                           ;-- counter for the colors array

;-- draw the color filled polygons
  do i=0, dimsizes(DEU_counties)-1  
     do ll=0, dimsizes(wcounties)-1
        if(names3(i) .eq. state_counties(ll)) then
           print("Draw county outline:  "+names3(i))
           startSegment = geometry(i, geom_segIndex)
           numSegments  = geometry(i, geom_numSegs)
           do seg=startSegment, startSegment+numSegments-1
              startPT = segments(seg, segs_xyzIndex)
              endPT   = startPT + segments(seg, segs_numPnts) - 1
              gnres@gsLineColor = "black"
              gnres@tfPolyDrawOrder = "PostDraw"
              poly(npl)  = gsn_add_polyline(wks, plot, lon(startPT:endPT),lat(startPT:endPT), gnres)
              npl = npl + 1
           end do
           j=j+1
        end if
     end do
  end do

  return(poly)            ;-- return the polygon array
end

;----------------------------------------------------------------------
;-- Function:  avg_by_county(...)
;--                ->  compute the average of the data for each county
;--------------------------------------------------------------------
undef("avg_by_county")
;--------------------------------------------------------------------
function avg_by_county(wks,plot,data,f,state_name,wcounties,levels,colors)
;--------------------------------------------------------------------
local f, segments, geometry, segsDims, geomDims, geom_segIndex, geom_numSegs, segs_xyzIndex, \
      segs_numPnts, numFeatures, i, lat, lon, startSegment, numSegments, seg, \
      startPT, endPT, dims, minlat, maxlat, minlon, maxlon
begin
  getvalues plot
      "mpLeftCornerLatF"      :   minlat                        ;-- retrieve map min lat
      "mpRightCornerLatF"     :   maxlat                        ;-- retrieve map max lat
      "mpLeftCornerLonF"      :   minlon                        ;-- retrieve map min lon
      "mpRightCornerLonF"     :   maxlon                        ;-- retrieve map max lon
  end getvalues

;-- read shapefile data
  geomDims      =  getfilevardimsizes(f,"geometry")
  numFeatures   =  geomDims(0)
  segments      =  f->segments
  geometry      =  f->geometry
  segsDims      =  dimsizes(segments)
  geom_segIndex =  f@geom_segIndex
  geom_numSegs  =  f@geom_numSegs
  segs_xyzIndex =  f@segs_xyzIndex
  segs_numPnts  =  f@segs_numPnts
  lat           =  f->y
  lon           =  f->x

  dims          =  dimsizes(data)
  nlat          =  dims(0)
  nlon          =  dims(1)
  lat1d         =  ndtooned(conform_dims((/nlat,nlon/),data&lat,0))
  lon1d         =  ndtooned(conform_dims((/nlat,nlon/),data&lon,1))
  nlatlon       =  dimsizes(lat1d)
  ii_latlon     =  ind(lat1d.ge.minlat.and.lat1d.le.maxlat.and.lon1d.ge.minlon.and.lon1d.le.maxlon)
  nii_latlon    =  dimsizes(ii_latlon)

;-- grab the indexes containing the counties
  states        =  f->NAME_1                                ;-- state names reference
  names3        =  f->NAME_3                                ;-- county names copied from DEU_adm3.shp
  DEU_counties  =  ind(names3.ne."")                           ;-- read all county names
 
;-- get state_name counties
  wc=new(dimsizes(names3),typeof(names3))
  if(.not.isatt(wc,"_FillValue")) then
    wc@_FillValue = default_fillvalue(typeof(names3))          ;-- make sure "wc" has a missing value
  end if
  n=0
  do m=0,dimsizes(names3)-1
    if(states(m).eq.state_name) then
       wc(n) = names3(m)                                       ;-- get all relevant counties
    else
       wc(n) = default_fillvalue(typeof(names3))               ;-- set all other to missing value
    end if
    n=n+1
  end do
  wcounties = ind(.not. ismissing(wc))
  state_counties = new(dimsizes(wcounties),string)
  do jj=0,dimsizes(wcounties)-1
     state_counties(jj) = names3(wcounties(jj))                ;-- get the names of the relevant counties
  end do

;-- create array to hold new data mask and averaged data
  data_mask_1d  = new(nlatlon,integer)
  if(.not.isatt(data,"_FillValue")) then
    data@_FillValue = default_fillvalue(typeof(data))          ;-- make sure "data" has a missing value
  end if
  data_1d  = ndtooned(data)                                    ;-- convert data to 1D array
  data_avg = new(dimsizes(DEU_counties),typeof(data),data@_FillValue)

  gnres    = True                                              ;-- polygon resource list
  nfill    = dimsizes(colors)

  do i=0,dimsizes(DEU_counties)-1
     do ll=0,dimsizes(wcounties)-1
        if (names3(i) .eq. state_counties(ll)) then     
           data_mask_1d = 0                                          ; Be sure to reset to 0 for every county
           startSegment = geometry(DEU_counties(i), geom_segIndex)
           numSegments  = geometry(DEU_counties(i), geom_numSegs)
           do seg=startSegment, startSegment+numSegments-1
              startPT = segments(seg, segs_xyzIndex)
              endPT   = startPT + segments(seg, segs_numPnts) - 1
              do n=0,nii_latlon-1
                 nn = ii_latlon(n)                                   ; Get index of point we're checking
                 if(lat1d(nn).lt.min(lat(startPT:endPT)).or.lat1d(nn).gt.max(lat(startPT:endPT)).or.\
                    lon1d(nn).lt.min(lon(startPT:endPT)).or.lon1d(nn).gt.max(lon(startPT:endPT))) 
                    continue
                 end if
                 if(gc_inout(lat1d(nn),lon1d(nn),lat(startPT:endPT),lon(startPT:endPT))) then
                    data_mask_1d(nn) = 1    ; This point is inside this county
                 end if
              end do
           end do

           ndm = num(data_mask_1d.eq.1)

           ;-- calculate the averages
           if(ndm.gt.0) then
             data_avg(i)  = avg(where(data_mask_1d.eq.1,data_1d,data_1d@_FillValue))
             print("-----------------------------------------------------------------")
             print((ll+1)+": Inspecting "+state_name+" county '" + names3(DEU_counties(i)) + "'...")
             print("     "+ndm + " data values found --> average = " + data_avg(i))
           end if
        end if
     end do
  end do

  return(data_avg)                               ;-- return data averages for each county
end

;----------------------------------------------------------------------
;-- Function:  avg_by_state(...)
;--                ->  compute the average of the data for each state
;--------------------------------------------------------------------
undef("avg_by_state")
;--------------------------------------------------------------------
function avg_by_state(wks,plot,data,shpf,state_names,levels,colors)
;--------------------------------------------------------------------
local segments, geometry, segsDims, geomDims, geom_segIndex, geom_numSegs, segs_xyzIndex, \
      segs_numPnts, numFeatures, i, lat, lon, startSegment, numSegments, seg, \
      startPT, endPT, dims, minlat, maxlat, minlon, maxlon
begin

  getvalues plot
      "mpLeftCornerLatF"      :   minlat
      "mpRightCornerLatF"     :   maxlat
      "mpLeftCornerLonF"      :   minlon
      "mpRightCornerLonF"     :   maxlon
  end getvalues
  
;-- read shapefile data
  geomDims      = getfilevardimsizes(shpf,"geometry")
  numFeatures   = geomDims(0)
  segments      = shpf->segments
  geometry      = shpf->geometry
  segsDims      = dimsizes(segments)
  geom_segIndex = shpf@geom_segIndex
  geom_numSegs  = shpf@geom_numSegs
  segs_xyzIndex = shpf@segs_xyzIndex
  segs_numPnts  = shpf@segs_numPnts
  lat           = shpf->y
  lon           = shpf->x

  dims          =  dimsizes(data)
  nlat          =  dims(0)
  nlon          =  dims(1)
  lat1d         =  ndtooned(conform_dims((/nlat,nlon/),data&lat,0))
  lon1d         =  ndtooned(conform_dims((/nlat,nlon/),data&lon,1))
  nlatlon       =  dimsizes(lat1d)
  ii_latlon     =  ind(lat1d.ge.minlat.and.lat1d.le.maxlat.and.lon1d.ge.minlon.and.lon1d.le.maxlon)
  nii_latlon    =  dimsizes(ii_latlon)

;-- grab the indexes containing the states
  names1        = state_names                             ;--  state names copied from DEU_adm1.shp
  ga_states     = ind(names1.ne."")

;-- create array to hold new data mask and averaged data
  data_mask_1d  = new(nlatlon,integer)
  
  if(.not. isatt(data,"_FillValue")) then
     data@_FillValue = default_fillvalue(typeof(data))    ;-- make sure "data" has a missing value
  end if
  
  data_1d       = ndtooned(data)                          ;-- this is for the averages computation later
  data_avg      = new(dimsizes(ga_states),typeof(data),data@_FillValue)  
                                                          ;-- array to hold data averages for each state
  gnres         = True                                    ;-- polygon resource list
  nfill         = dimsizes(colors)

  do i=0, dimsizes(ga_states)-1
     data_mask_1d = 0                                     ;-- be sure to reset to 0 for every state
     startSegment = geometry(ga_states(i), geom_segIndex)
     numSegments  = geometry(ga_states(i), geom_numSegs)  ;-- some states have multiple segments
     do seg=startSegment, startSegment+numSegments-1
        startPT = segments(seg, segs_xyzIndex)
        endPT   = startPT + segments(seg, segs_numPnts) - 1
        do n=0,nii_latlon-1                         ;-- loop through each point if it's in this state
           nn = ii_latlon(n)                        ;-- get index of point we're checking
           ;-- don't check this point if it doesn't fall within this state
           if(lat1d(nn).lt.min(lat(startPT:endPT)).or.lat1d(nn).gt.max(lat(startPT:endPT)).or.\
              lon1d(nn).lt.min(lon(startPT:endPT)).or.lon1d(nn).gt.max(lon(startPT:endPT))) 
              continue
           end if
           if(gc_inout(lat1d(nn),lon1d(nn),lat(startPT:endPT),lon(startPT:endPT))) then
              data_mask_1d(nn) = 1                   ;-- this point is inside this state
           end if
        end do
    end do

    ndm = num(data_mask_1d.eq.1)                     ;-- count number of points found in this state

;-- compute the averages
    if(ndm.gt.0) then
       data_avg(i)  =  avg(where(data_mask_1d.eq.1,data_1d,data_1d@_FillValue))
       print("-----------------------------------------------------------------")
       print((i+1)+": Inspecting Germany state '" + names1(ga_states(i)) + "'...")
    end if
    
    print("     "+ndm + " data values found --> average = " + data_avg(i))
       
  end do

  return(data_avg)                                    ;-- return data averages for each state
end

;----------------------------------------------------------------------
;-- MAIN
;----------------------------------------------------------------------
begin
  start_date = toint(systemfunc("date +%s"))
  state_name = (/"Schleswig-Holstein","Hamburg","Bremen", "Niedersachsen","Mecklenburg-Vorpommern"/)
  print("")
  print("Selected state of Germany:  "+state_name)
  print("")

  mminlat  =  51.0
  mmaxlat  =  55.5
  mminlon  =  6.5
  mmaxlon  =  14.75

  country_border  =  True                                 ;-- default: draw country border line
  states_border   =  True                                 ;-- default: draw states border lines
  counties_border =  True                                 ;-- default: draw counties border lines
  subregion       =  True
  
;-- shapefile containing Germany states and counties
  shapefile_dir  = "$HOME/data/Shapefiles/DEU_adm/"       ;-- directory containing the shapefiles
  shp_name       = "DEU_adm3.shp"                         ;-- shapefile to be used
  shp_fname      = shapefile_dir+shp_name                 ;-- path of shapefile
  shpf3          = addfile(shp_fname,"r")                 ;-- open shapefile
  county_names   = shpf3->NAME_3                          ;-- county names
  states         = shpf3->NAME_1                          ;-- state names
  shplon         = shpf3->x                               ;-- longitudes
  shplat         = shpf3->y                               ;-- latitudes

;-- Germany borderline coordinates
  DEU_minlat     =  min(shplat)-0.1
  DEU_maxlat     =  max(shplat)+0.1
  DEU_minlon     =  min(shplon)-0.1
  DEU_maxlon     =  max(shplon)+0.1

;-- generate dummy data (we need higher resolution for regional sections)
  nlat           =  150                                    ;-- number of lat values
  nlon           =  150                                    ;-- number of lon values
  lat            =  fspan(DEU_minlat,DEU_maxlat,nlat)      ;-- generate lat dimension data
  lon            =  fspan(DEU_minlon,DEU_maxlon,nlon)      ;-- generate lon dimension data
  lat@units      = "degrees_north"                         ;-- lat dimension units attribute
  lon@units      = "degrees_east"                          ;-- lon dimension units attribute
  var            =  generate_2d_array(25,25,-15,20,100,(/nlat,nlon/))  ;-- generate dummy data
  var!0          = "lat"                                   ;-- data lat dimension name
  var!1          = "lon"                                   ;-- data lon dimension name
  var&lat        =  lat                                    ;-- lat dimension data
  var&lon        =  lon                                    ;-- lon dimension data
  var@units      = "C"                                     ;-- data units
  var@_FillValue = -9999.9
  
  delta_x = (DEU_maxlon-DEU_minlon)/nlon
  delta_y = (DEU_maxlat-DEU_minlat)/nlat
  
;-- open a workstation
  wks = gsn_open_wks("png","plot_DEU_adm3_avg_over_counties_COAST")

;-- set resources
  res                       =  True
  res@gsnDraw               =  False                 ;-- don't draw plot yet
  res@gsnFrame              =  False                 ;-- don't advance frame
  res@gsnAddCyclic          =  False                 ;-- no cyclic point
  res@gsnRightString        = ""
  
  res@mpProjection          = "Mercator"             ;-- use Mercator projection
  res@mpLimitMode           = "Corners"              ;-- map limit mode
  if(isvar("subregion")) then                        ;-- is 'subregion' on command line?
     res@mpLeftCornerLatF     =  mminlat             ;-- min lat
     res@mpRightCornerLatF    =  mmaxlat             ;-- max lat
     res@mpLeftCornerLonF     =  mminlon             ;-- min lon
     res@mpRightCornerLonF    =  mmaxlon             ;-- max lon
  else
     res@mpLeftCornerLatF     =  DEU_minlat          ;-- min lat
     res@mpRightCornerLatF    =  DEU_maxlat          ;-- max lat
     res@mpLeftCornerLonF     =  DEU_minlon          ;-- min lon
     res@mpRightCornerLonF    =  DEU_maxlon          ;-- max lon
  end if
  res@mpDefaultFillColor    =  16                    ;-- draw land in gray
  res@mpOutlineOn           =  True                  ;-- draw map outlines
  res@mpDataBaseVersion     = "HighRes"              ;-- use HighRes map
  res@mpDataResolution      = "Fine"                 ;-- we need a finer resolution

  res@tiMainFontHeightF     =  0.018                 ;-- title font size

  res@pmTickMarkDisplayMode = "Always"               ;-- better tickmark labels

  res@vpHeightF             =  0.72
  res@vpWidthF              =  1.0
  res@vpXF                  =  0.01
  res@vpYF                  =  0.84

  res@tmYLLabelFontHeightF  =  0.013
  res@tmXBLabelFontHeightF  =  0.013 
  res@tmXBMajorLengthF      =  0.01
  res@tmYLMajorLengthF      =  0.01
  
  res@tiMainString         = "  "

  res@cnFillOn              =  True                  ;-- contour fill on
  res@cnLinesOn             =  False                 ;-- turn off contour lines
  res@cnLineLabelsOn        =  False                 ;-- turn off contour line labels
  res@cnLevelSelectionMode  = "ManualLevels"         ;-- set levels
  res@cnMinLevelValF        =  min(var)              ;-- min values
  res@cnMaxLevelValF        =  max(var)              ;-- max values
  res@cnLevelSpacingF       =  0.5                   ;-- increment value
  res@cnFillPalette         = "rainbow"              ;-- colormap
  res@cnMissingValFillColor =  -1                    ;-- set missing fill color to 100% transparency

  res@lbLabelBarOn          =  False
  
  plot_orig = gsn_csm_contour_map(wks,var,res)       ;-- create contour plot to retrieve the
                                                     ;-- levels and colors values, but don't draw it  

;-- this gives us the colors and levels to use for the filled polygons
  getvalues plot_orig@contour
     "cnLevels"               : levels                ;-- retrieve levels
     "cnFillColors"           : colors                ;-- retrieve colors
     "cnInfoLabelFontHeightF" : font_h                ;-- retrieve font height
  end getvalues

  map = setColorContourClear(plot_orig,min(var),max(var))
  
;-----------------------------------
;-- 1. State - Schleswig-Holstein
;-----------------------------------
     wc1 = new(dimsizes(county_names),typeof(county_names))     ;-- assign array for the selected counties
     if(.not.isatt(wc1,"_FillValue")) then
        wc1@_FillValue = default_fillvalue(typeof(county_names)) ;-- make sure "wc" has a missing value
     end if
     n=0
     do m=0,dimsizes(county_names)-1
        if(states(m).eq.state_name(0)) then
           wc1(n) = county_names(m)                             ;-- get counties of the state
        else
           wc1(n) = default_fillvalue(typeof(county_names))     ;-- set other counties to missing value
        end if
        n=n+1
     end do
     wcounties1 = ind(.not. ismissing(wc1))                     ;-- indices of counties
  
;-- calculate the averages
     var_avg1 = avg_by_county(wks, map, var, shpf3, state_name(0), wcounties1, levels, colors)

     print("--------------------------------------------------")
     print("Data values: " + num(.not. ismissing(var_avg1))+ "  Missing values: " + num(ismissing(var_avg1)))
  
;-- get the correct color indices for the averaged data
     col_avg1 = new(dimsizes(wc1),integer)                ;-- assign new color map
     do i=0,dimsizes(wc1)-1
        if(ismissing(var_avg1(i))) then
          col_avg1(i) = res@cnMissingValFillColor        ;-- if missing value use cnMissingValFillColor
       else
         do j=0,dimsizes(levels)-1
            if (var_avg1(i).lt.levels(0)) then
                col_avg1(i) = colors(0)                 ;-- values less than min(levels)
            else if(var_avg1(i).ge.levels(dimsizes(levels)-1)) then
                col_avg1(i) = colors(dimsizes(colors)-1) ;-- values greater than max(levels)
            else if(var_avg1(i).ge.levels(j).and.var_avg1(i).lt.levels(j+1)) then
                col_avg1(i) = colors(j+1)               ;-- values in between
            end if
            end if
            end if
         end do
      end if
    end do
    print("--------------------------------------------------")

;-- draw only the colored data averages in the counties polygons of selected state
    dum_poly1 = add_shapefile_polygons(wks, map, state_name(0), col_avg1, shpf3, True)
    print("added polygons to the plot ... done")
    
;-- draw all county polylines on top
    if(counties_border) then
       dum_polyl1 = add_shapefile_polylines(wks, map, state_name(0), col_avg1, shpf3, True)
       print("added polylines to the plot ... done")
    end if

;-----------------------------------
;-- 2. State - Hamburg
;-----------------------------------
     wc2 = new(dimsizes(county_names),typeof(county_names))     ;-- assign array for the selected counties
     if(.not.isatt(wc2,"_FillValue")) then
        wc2@_FillValue = default_fillvalue(typeof(county_names)) ;-- make sure "wc" has a missing value
     end if
     n=0
     do m=0,dimsizes(county_names)-1
        if(states(m).eq.state_name(1)) then
           wc2(n) = county_names(m)                             ;-- get counties of the state
        else
           wc2(n) = default_fillvalue(typeof(county_names))     ;-- set other counties to missing value
        end if
        n=n+1
     end do
     wcounties2 = ind(.not. ismissing(wc2))                     ;-- indices of counties
  
;-- calculate the averages
     var_avg2 = avg_by_county(wks, map, var, shpf3, state_name(1), wcounties2, levels, colors)

     print("--------------------------------------------------")
     print("Data values: " + num(.not. ismissing(var_avg2))+ "  Missing values: " + num(ismissing(var_avg2)))
  
;-- get the correct color indices for the averaged data
     col_avg2 = new(dimsizes(wc2),integer)                ;-- assign new color map
     do i=0,dimsizes(wc2)-1
        if(ismissing(var_avg2(i))) then
          col_avg2(i) = res@cnMissingValFillColor        ;-- if missing value use cnMissingValFillColor
       else
         do j=0,dimsizes(levels)-1
            if (var_avg2(i).lt.levels(0)) then
                col_avg2(i) = colors(0)                 ;-- values less than min(levels)
            else if(var_avg2(i).ge.levels(dimsizes(levels)-1)) then
                col_avg2(i) = colors(dimsizes(colors)-1) ;-- values greater than max(levels)
            else if(var_avg2(i).ge.levels(j).and.var_avg2(i).lt.levels(j+1)) then
                col_avg2(i) = colors(j+1)               ;-- values in between
            end if
            end if
            end if
         end do
      end if
    end do
    print("--------------------------------------------------")

;-- draw only the colored data averages in the counties polygons of selected state
    dum_poly2 = add_shapefile_polygons(wks, map, state_name(1), col_avg2, shpf3, True)
    print("added polygons to the plot ... done")
    
;-- draw all county polylines on top
    if(counties_border) then
       dum_polyl2 = add_shapefile_polylines(wks, map, state_name(1), col_avg2, shpf3, True)
       print("added polylines to the plot ... done")
    end if

;-----------------------------------
;-- 3. State - Bremen
;-----------------------------------
     wc3 = new(dimsizes(county_names),typeof(county_names))     ;-- assign array for the selected counties
     if(.not.isatt(wc3,"_FillValue")) then
        wc3@_FillValue = default_fillvalue(typeof(county_names)) ;-- make sure "wc" has a missing value
     end if
     n=0
     do m=0,dimsizes(county_names)-1
        if(states(m).eq.state_name(2)) then
           wc3(n) = county_names(m)                             ;-- get counties of the state
        else
           wc3(n) = default_fillvalue(typeof(county_names))     ;-- set other counties to missing value
        end if
        n=n+1
     end do
     wcounties3 = ind(.not. ismissing(wc3))                     ;-- indices of counties
  
;-- calculate the averages
     var_avg3 = avg_by_county(wks, map, var, shpf3, state_name(2), wcounties3, levels, colors)

     print("--------------------------------------------------")
     print("Data values: " + num(.not. ismissing(var_avg3))+ "  Missing values: " + num(ismissing(var_avg3)))
  
;-- get the correct color indices for the averaged data
     col_avg3 = new(dimsizes(wc3),integer)                ;-- assign new color map
     do i=0,dimsizes(wc3)-1
        if(ismissing(var_avg3(i))) then
          col_avg3(i) = res@cnMissingValFillColor        ;-- if missing value use cnMissingValFillColor
       else
         do j=0,dimsizes(levels)-1
            if (var_avg3(i).lt.levels(0)) then
                col_avg3(i) = colors(0)                 ;-- values less than min(levels)
            else if(var_avg3(i).ge.levels(dimsizes(levels)-1)) then
                col_avg3(i) = colors(dimsizes(colors)-1) ;-- values greater than max(levels)
            else if(var_avg3(i).ge.levels(j).and.var_avg3(i).lt.levels(j+1)) then
                col_avg3(i) = colors(j+1)               ;-- values in between
            end if
            end if
            end if
         end do
      end if
    end do
    print("--------------------------------------------------")

;-- draw only the colored data averages in the counties polygons of selected state
    dum_poly3 = add_shapefile_polygons(wks, map, state_name(2), col_avg3, shpf3, True)
    print("added polygons to the plot ... done")
    
;-- draw all county polylines on top
    if(counties_border) then
       dum_polyl3 = add_shapefile_polylines(wks, map, state_name(2), col_avg3, shpf3, True)
       print("added polylines to the plot ... done")
    end if

;-----------------------------------
;-- 4. State - Niedersachsen
;-----------------------------------
     wc4 = new(dimsizes(county_names),typeof(county_names))     ;-- assign array for the selected counties
     if(.not.isatt(wc4,"_FillValue")) then
        wc4@_FillValue = default_fillvalue(typeof(county_names)) ;-- make sure "wc" has a missing value
     end if
     n=0
     do m=0,dimsizes(county_names)-1
        if(states(m).eq.state_name(3)) then
           wc4(n) = county_names(m)                             ;-- get counties of the state
        else
           wc4(n) = default_fillvalue(typeof(county_names))     ;-- set other counties to missing value
        end if
        n=n+1
     end do
     wcounties4 = ind(.not. ismissing(wc4))                     ;-- indices of counties
  
;-- calculate the averages
     var_avg4 = avg_by_county(wks, map, var, shpf3, state_name(3), wcounties4, levels, colors)

     print("--------------------------------------------------")
     print("Data values: " + num(.not. ismissing(var_avg4))+ "  Missing values: " + num(ismissing(var_avg4)))
  
;-- get the correct color indices for the averaged data
     col_avg4 = new(dimsizes(wc4),integer)                ;-- assign new color map
     do i=0,dimsizes(wc4)-1
        if(ismissing(var_avg4(i))) then
          col_avg4(i) = res@cnMissingValFillColor        ;-- if missing value use cnMissingValFillColor
       else
         do j=0,dimsizes(levels)-1
            if (var_avg4(i).lt.levels(0)) then
                col_avg4(i) = colors(0)                 ;-- values less than min(levels)
            else if(var_avg4(i).ge.levels(dimsizes(levels)-1)) then
                col_avg4(i) = colors(dimsizes(colors)-1) ;-- values greater than max(levels)
            else if(var_avg4(i).ge.levels(j).and.var_avg4(i).lt.levels(j+1)) then
                col_avg4(i) = colors(j+1)               ;-- values in between
            end if
            end if
            end if
         end do
      end if
    end do
    print("--------------------------------------------------")

;-- draw only the colored data averages in the counties polygons of selected state
    dum_poly4 = add_shapefile_polygons(wks, map, state_name(3), col_avg4, shpf3, True)
    print("added polygons to the plot ... done")
    
;-- draw all county polylines on top
    if(counties_border) then
       dum_polyl4 = add_shapefile_polylines(wks, map, state_name(3), col_avg4, shpf3, True)
       print("added polylines to the plot ... done")
    end if

;-----------------------------------
;-- 5. State - Mecklenburg-Vorpommern
;-----------------------------------
     wc5 = new(dimsizes(county_names),typeof(county_names))     ;-- assign array for the selected counties
     if(.not.isatt(wc5,"_FillValue")) then
        wc5@_FillValue = default_fillvalue(typeof(county_names)) ;-- make sure "wc" has a missing value
     end if
     n=0
     do m=0,dimsizes(county_names)-1
        if(states(m).eq.state_name(4)) then
           wc5(n) = county_names(m)                             ;-- get counties of the state
        else
           wc5(n) = default_fillvalue(typeof(county_names))     ;-- set other counties to missing value
        end if
        n=n+1
     end do
     wcounties5 = ind(.not. ismissing(wc5))                     ;-- indices of counties
  
;-- calculate the averages
     var_avg5 = avg_by_county(wks, map, var, shpf3, state_name(4), wcounties5, levels, colors)

     print("--------------------------------------------------")
     print("Data values: " + num(.not. ismissing(var_avg5))+ "  Missing values: " + num(ismissing(var_avg5)))
  
;-- get the correct color indices for the averaged data
     col_avg5 = new(dimsizes(wc5),integer)                ;-- assign new color map
     do i=0,dimsizes(wc5)-1
        if(ismissing(var_avg5(i))) then
          col_avg5(i) = res@cnMissingValFillColor        ;-- if missing value use cnMissingValFillColor
       else
         do j=0,dimsizes(levels)-1
            if (var_avg5(i).lt.levels(0)) then
                col_avg5(i) = colors(0)                 ;-- values less than min(levels)
            else if(var_avg5(i).ge.levels(dimsizes(levels)-1)) then
                col_avg5(i) = colors(dimsizes(colors)-1) ;-- values greater than max(levels)
            else if(var_avg5(i).ge.levels(j).and.var_avg5(i).lt.levels(j+1)) then
                col_avg5(i) = colors(j+1)               ;-- values in between
            end if
            end if
            end if
         end do
      end if
    end do
    print("--------------------------------------------------")

;-- draw only the colored data averages in the counties polygons of selected state
    dum_poly5 = add_shapefile_polygons(wks, map, state_name(4), col_avg5, shpf3, True)
    print("added polygons to the plot ... done")
    
;-- draw all county polylines on top
    if(counties_border) then
       dum_polyl5 = add_shapefile_polylines(wks, map, state_name(4), col_avg5, shpf3, True)
       print("added polylines to the plot ... done")
    end if

;------------------------------------------
;-- draw the states borderlines of Germany
;------------------------------------------
  if(states_border) then
    shp_name1  = "DEU_adm1.shp"                      ;-- shapefile to be used
    shp_fname1 = shapefile_dir+shp_name1             ;-- path of shapefile
    sborder = gsn_add_shapefile_polylines(wks, map, shp_fname1, True)
    print("added states polylines to the plot ... done")
  end if
  
;------------------------------------------
;-- draw the borderline of Germany
;------------------------------------------
  if(country_border) then
    shp_name0  = "DEU_adm0.shp"                      ;-- shapefile to be used
    shp_fname0 = shapefile_dir+shp_name0             ;-- path of shapefile
    cborder = gsn_add_shapefile_polylines(wks, map, shp_fname0, True)
    print("added country polylines to the plot ... done")
  end if    

;------------------------------------------
;-- add a common labelbar
;------------------------------------------
  lbres                      =  True
  lbres@lbPerimOn            =  False               ;-- don't draw labelbar boxes
  lbres@lbOrientation        = "Horizontal"         ;-- labelbar orientation
  lbres@vpWidthF             =  0.7                 ;-- width of the labelbar
  lbres@vpHeightF            =  0.08                ;-- height of the labelbar
  lbres@lbLabelFontHeightF   =  0.012               ;-- labelbar label font height
  lbres@lbLabelAlignment     = "InteriorEdges"      ;-- labelbar label alignment
  lbres@lbMonoFillPattern    =  True                ;-- labelbar solid fill
  lbres@lbFillColors         =  colors              ;-- labelbar colors (must be RGB triplets)
  labels  = "" + levels                             ;-- set labels
  nlevels = dimsizes(levels)                        ;-- number of levels
  gsn_labelbar_ndc (wks,nlevels+1,labels,0.16,0.084,lbres) ;-- add labelbar
  
;------------------------------------------
;-- add title strings
;------------------------------------------
  title0 = "Germany"
  title1 = "data averaged over the counties"
  title2 = "(grid:  dlon="+sprintf("%5.3f",delta_x)+"~S~o~N~  dlat="+sprintf("%5.3f",delta_y)+"~S~o~N~)"
  
  names = state_name
  if(dimsizes(state_name).gt.1) then
     names := "    "+state_name(0)
     do mm=1,dimsizes(state_name)-1
        if(dimsizes(state_name).gt.3.and.mm.eq.3) then
           names := names + ",~C~" + state_name(mm)
        else
           names := names + ", " + state_name(mm)
        end if
     end do
  end if

  tires                   =  True                   ;-- text resources title string
  tires@txFontHeightF     =  0.024                  ;-- text font size
  res@txFontThicknessF    =  1.8                    ;-- bold
  tires@txJust            = "CenterCenter"          ;-- text justification
  gsn_text_ndc(wks,names,0.5, 0.950, tires)         ;-- center upper title string
  tires@txFontHeightF     =  0.014                  ;-- text font size
  gsn_text_ndc(wks,title0, 0.5, 0.905, tires)       ;-- center middle title string
  gsn_text_ndc(wks,title1, 0.5, 0.885, tires)       ;-- center middle title string
  tires@txFontHeightF     =  0.012                  ;-- text font size
  gsn_text_ndc(wks,title2, 0.5, 0.865, tires)       ;-- center lower title string
;------------------------------------------
;-- add units to labelbar and the copyright string
;------------------------------------------
  tires@txJust            = "BottomRight"           ;-- text justification
  tires@txFontHeightF     =  0.012                  ;-- make font size smaller
  gsn_text_ndc(wks,"Temperature  [~S~o~N~C]", 0.6, 0.005, tires) ;-- plot units string
  gsn_text_ndc(wks,"~F35~c ~F21~~N~DKRZ",    0.92, 0.005, tires) ;-- plot copyright info 
;------------------------------------------
;-- create the plot and advance the frame
;------------------------------------------
  draw(map)
  frame(wks)
;------------------------------------------ 
  end_date = toint(systemfunc("date +%s"))
  print("Elapsed time:  "+(end_date-start_date)+"s")
end

Artikelaktionen