Re: Calculation and plotting of the climate zones after the Köppen & Geiger climate classification

Hi Stephan,

the do loops with loc_lat and loc_lon are wro= ng. You have to change it to

<= tt>do loc_lat = 0,dimsizes(lat2d(:,0))-1
do loc_l= on = 0,dimsizes(lon2d(0,:))-1

To define the var= iable Klimaklasse with all attributes and dimensions of
the variable= Niederschlag you shouldn't use 'new' but copy it:

Klimaklasse = Niederschlag
The changes can be found in your attached script.

Bye, Karin

Am 21.06.14 19:15, schrieb Stephan Herrmann:

Hello everybody,

 

I want to calculate an= d plot the climate zones after the Köppen & Geiger climate classif= ication with the following NCL code:

 

load "$NCARG_ROOT/lib/= ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscr= ipts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/cont= ributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl= "

begin

f = addfile("MPIM-REMO_EH5OMA1Br3_HR_04= 4_BOL_1961-1990_167_minus_1_hour_Celsiustt_ymonmean.nc","r")
Temperat= ur = f->var167(:,:,:,:)

latsize = getfilevardimsizes(f,= "rlat")
lonsize = getfilevardimsizes(f,"rlon")
rotpole = f-= >rotated_pole        ; need attribute= s for correct remap

ll = asciiread("ll.out",(/latsize,lonsiz= e,2/),"float")
lat2d = ll(::,:,0)
lon2d = ll(::,:,1)
=
Temperatur@lat2d = lat2d
Temperatur@lon2d = lon2d
<= br /> f = addfile("MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minu= s_1_hour_ymonsum_mean_30_years.nc","r")
Niederschlag = f->var142= (:,:,:)

latsize = getfilevardimsizes(f,"rlat")
lonsize= = getfilevardimsizes(f,"rlon")
rotpole = f->rotated_pole = ;       ; need attributes for correct remap
ll = asciiread("ll.out",(/latsize,lonsize,2/),"float")
= lat2d = ll(::,:,0)
lon2d = ll(::,:,1)

Niederschlag@= lat2d = lat2d
Niederschlag@lon2d = lon2d

Klimaklasse= = new(dimsizes(Niederschlag),float)
    
Klim= aklasse = 99999.9

do loc_lat = 0,max(dimsizes(lat2d))-1 do loc_lon = 0,max(dimsizes(lon2d))-1

  temp = n= ew(12,typeof(Temperatur))
  temp = Temperatur(:,0,loc_lat,loc_= lon)

  nied = new(12,typeof(Niederschlag))
 = nied = Niederschlag(:,loc_lat,loc_lon)

  summe_nied = = 0.0
  mean_temp = 0.0
  do i = 0,11
&nb= sp;   summe_nied = summe_nied + nied(i)
   = ; mean_temp = mean_temp + temp(i)
  end do
  mean_t= emp = mean_temp / 12.0

;**** Klimaklassifikation nach Koeppe= n-Geiger

 hk = 6

 sommer_nied = 0.0=

 do k = min((/mod(3 + hk, 12), mod(8 + hk, 12)/)), max= ((/mod(3 + hk, 12), mod(8 + hk, 12)/)) ; oder hydrologisch Mai bis Oktober?=
    sommer_nied = sommer_nied + nied(k)
 = ;end do

 if(summe_nied * 0.7 .lt. sommer_nied) then
=     regen = 14.0       = ;       ; Sommerregen
 else if(sum= me_nied * 0.3 .gt. sommer_nied) then
    regen = 0.0=             &nb= sp;  ; Winterregen
 else
    regen = = 7.0           &n= bsp;   ; gleichmäßige Niederschlagsverteilung
&n= bsp;end if
 end if

 warm = 0.0
  =
 do k=0,11
    if(temp(k) .gt. 10.0) the= n
       warm = warm + 1.0
&nbs= p;   end if
 end do

 if((min(temp(:)= ) .gt. 18.0) .and. (summe_nied / 10.0 .gt. 2.0 * (mean_temp + regen))) then=
    klasse = "A"      = ;         ; Tropisches Regenwald- o= der Savannenklima ohne Winter

    if(min(nied(:= )) .gt. 60.0) then
       klasse = kl= asse + "f"   ; Tropisches Regenwaldklima
   = else if (min(nied(:)) .gt. 0.04 * (2400.0 - summe_nied)) then ;aendern sie= he Zettel
       klasse = klasse + "m= "   ; Mittelform, Regenwaldklima trotz einer Trockenzeit
&n= bsp;   else
       klasse == klasse + "w"   ; Savannenklima
    end if     end if

 else if(summe_nied / 10.0= .lt. 2.0 * (mean_temp + regen)) then
   klasse = "B"&nbs= p;            &= nbsp; ; Trockenklima

   if(summe_nied / 10.0 .gt. me= an_temp + regen) then
        klas= se = klasse + "S"   ; Steppenklima
   else
        klasse = klasse + "W" &= nbsp; ; Wuestenklima
    end if
   &nb= sp;
    if(mean_temp .gt. 18.0) then
  = ;     klasse = klasse + "h"   ; heiß<= br />     else if(mean_temp .le. 18.0) then
 &nbs= p;     klasse = klasse + "k"   ; kalt
=     end if
  end if

  else if(= (min(temp(:)) .gt. -3.0) .and. (min(temp(:)) .lt. 18.0) .and. (summe_nied /= 10.0 .gt. 2.0 * (mean_temp + regen))) then
    klasse= = "C"           &= nbsp;   ; Warm-gemaeßigtes Klima

  &= nbsp; if(min((/nied(mod(11 + hk, 12)),nied(mod(0 + hk, 12)),nied(mod(1 + hk= , 12))/)) .lt. 0.1 * max((/nied(mod(5 + hk, 12)),nied(mod(6 + hk, 12)),nied= (mod(7 + hk, 12))/))) then
      

klasse = klasse + "w= "   ; warmes wintertrockenes Klima

 

    els= e if((min((/nied(mod(5 + hk, 12)),nied(mod(6 + hk, 12)),nied(mod(7 + hk, 12= ))/)) .lt. 40.0) .and. (max((/nied(mod(11 + hk, 12)),nied(mod(0 + hk, 12)),= nied(mod(1 + hk, 12))/) / 3.0) .gt. min((/nied(mod(5 + hk, 12)),nied(mod(6 = + hk, 12)),nied(mod(7 + hk, 12))/)))) then
    &n= bsp;  klasse = klasse + "s"   ; warmes sommertrockenes Kli= ma
    else
      = klasse = klasse + "f"   ; feuchtgemaessigtes Klima
 = ;   end if
    end if

 &nb= sp;  if(max(temp(:)) .gt. 22.0) then
    &nb= sp;  klasse = klasse + "a"   ;heisse Sommer
 &n= bsp;  else if((max(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) then ;ab= hier neu einfuegen
       klasse = k= lasse + "b"   ; warme Sommer
    else if((ma= x(temp(:)) .lt. 22.0) .and. (warm .le. 3.0) .and. (warm .ge. 1.0)) then
       klasse = klasse + "c"  = ; kuehle Sommer
    end if
    e= nd if
    end if

  else if((min(temp= (:)) .lt. -3.0) .and. (max(temp(:)) .gt. 10.0)) then
  &nbs= p; klasse = "D"         &nbs= p;     ; Boreales oder Schnee-Wald-Klima
 &n= bsp;  
    if(min((/nied(mod(11 + hk, 12)),nied(m= od(0 + hk, 12)),nied(mod(1 + hk, 12))/)) .lt. 0.1 * max((/nied(mod(5 + hk, = 12)),nied(mod(6 + hk, 12)),nied(mod(7 + hk, 12))/))) then
  = ;     klasse = klasse + "w"   ; wintertrock= enkaltes Klima
    else
    =    klasse = klasse + "f"   ; winterfeuchtkaltes Klima=
    end if
    
 &nbs= p;  if(max(temp(:)) .gt. 22.0) then
    &nbs= p;  klasse = klasse + "a"   ; heiße Sommer
&n= bsp;   else if((max(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) th= en
       klasse = klasse + "b" =   ; warme Sommer
    else if((max(temp(:)) .lt. 2= 2.0) .and. (warm .le. 3.0) .and. (warm .ge. 1.0)) then
  &n= bsp;    klasse = klasse + "c"   ; kuehle Sommer<= br />     else if(min(temp(:)) .lt. -38.0) then
 =       klasse = klasse + "d"   ; streng= e Winter
    end if
    end if     end if
    end if
 &= nbsp;  
  else if(max(temp(:)) .lt. 10.0) then
 =    klasse = "E"        =        ; Schneeklima
    = ;
    if(max(temp(:)) .lt. 0.0) then
  = ;     klasse = klasse + "F"   ; Tundrenklim= a
    else
       = klasse = klasse + "T"   ; Klima ewigen Frostes
 &nbs= p;  end if
    
  end if
  = end if
  end if
  end if
  end if
  if (klasse .eq. "Af") then
     Klima= klasse(0,loc_lat,loc_lon) = 30.1
  else if (klasse .eq. "Am") = then
     Klimaklasse(0,loc_lat,loc_lon) = 29.1=
  else if (klasse .eq. "As") then
    = ; Klimaklasse(0,loc_lat,loc_lon) = 28.1
  else if (klasse .eq.= "Aw") then
     Klimaklasse(0,loc_lat,loc_lon) = = 27.1
  else if (klasse .eq. "BWk") then
  &n= bsp;  Klimaklasse(0,loc_lat,loc_lon) = 26.1
  else if (kl= asse .eq. "BWh") then
     Klimaklasse(0,loc_lat,= loc_lon) = 25.1
  else if (klasse .eq. "BSk") then
 = ;    Klimaklasse(0,loc_lat,loc_lon) = 24.1
  e= lse if (klasse .eq. "BSh") then
     Klimaklasse(= 0,loc_lat,loc_lon) = 23.1
  else if (klasse .eq. "Cfa") then      Klimaklasse(0,loc_lat,loc_lon) = 22.1
=   else if (klasse .eq. "Cfb") then
     Kli= maklasse(0,loc_lat,loc_lon) = 21.1
  else if (klasse .eq. "Cfc= ") then
     Klimaklasse(0,loc_lat,loc_lon) = 2= 0.1
  else if (klasse .eq. "Csa") then
   &= nbsp; Klimaklasse(0,loc_lat,loc_lon) = 19.1
  else if (klasse = .eq. "Csb") then
     Klimaklasse(0,loc_lat,loc_l= on) = 18.1
  else if (klasse .eq. "Csc") then
 &nbs= p;   Klimaklasse(0,loc_lat,loc_lon) = 17.1
  else if= (klasse .eq. "Cwa") then
     Klimaklasse(0,loc_= lat,loc_lon) = 16.1
  else if (klasse .eq. "Cwb") then
&= nbsp;    Klimaklasse(0,loc_lat,loc_lon) = 15.1
 = ; else if (klasse .eq. "Cwc") then
     Klimaklas= se(0,loc_lat,loc_lon) = 14.1
  else if (klasse .eq. "Dfa") the= n
     Klimaklasse(0,loc_lat,loc_lon) = 13.1   else if (klasse .eq. "Dfb") then
     = Klimaklasse(0,loc_lat,loc_lon) = 12.1
  else if (klasse .eq. "= Dfc") then
     Klimaklasse(0,loc_lat,loc_lon) = = 11.1
  else if (klasse .eq. "Dfd") then
  &n= bsp;  Klimaklasse(0,loc_lat,loc_lon) = 10.1
  else if (kl= asse .eq. "Dsa") then
     Klimaklasse(0,loc_lat,= loc_lon) = 9.1
  else if (klasse .eq. "Dsb") then
 =     Klimaklasse(0,loc_lat,loc_lon) = 8.1
  else= if (klasse .eq. "Dsc") then
     Klimaklasse(0,l= oc_lat,loc_lon) = 7.1
   else if (klasse .eq. "Dsd") then=
     Klimaklasse(0,loc_lat,loc_lon) = 6.1
  else if (klasse .eq. "Dwa") then
     Kl= imaklasse(0,loc_lat,loc_lon) = 5.1
  else if (klasse .eq. "Dwb= ") then
     Klimaklasse(0,loc_lat,loc_lon) = 4= .1
  else if (klasse .eq. "Dwc") then
   &n= bsp; Klimaklasse(0,loc_lat,loc_lon) = 3.1
  else if (klasse .e= q. "Dwd") then
     Klimaklasse(0,loc_lat,loc_lon= ) = 2.1
  else if (klasse .eq. "EF") then
  &n= bsp;  Klimaklasse(0,loc_lat,loc_lon) = 1.1
  else if (kla= sse .eq. "ET") then
     Klimaklasse(0,loc_lat,lo= c_lon) = 0.1
  end if
  end if
  end if<= br />   end if
  end if
  end if
  en= d if
  end if
  end if
  end if
&nbs= p; end if
  end if
  end if
  end if
=   end if
  end if
  end if
  end if<= br />   end if
  end if
  end if
  en= d if
  end if
  end if
  end if
&nbs= p; end if
  end if
  end if
  end if
=   end if
  end if

 
end do
= end do

wks = gsn_open_wks("ps","Koeppen")   &= nbsp;           ; open a = ps file
;plot = new(4,graphic)      &= nbsp;           &nbs= p;       ; create a plot array

g= sn_define_colormap(wks,"cb_9step")

res = True
res@gsnM= aximize = True
res@gsnAddCyclic = False
res@mpDataBaseVersi= on = "mediumres"
res@mpOutlineBoundarySets = "national"
&nb= sp;       
res@mpCenterLonF  = = rotpole@grid_north_pole_longitude
res@mpCenterLatF   = = rotpole@grid_north_pole_latitude -90
res@cnFillOn = True
= res@cnLinesOn = False
res@cnLineLabelsOn = False

res= @cnExplicitLegendLabelsOn = True
;res@lbLabelAlignment = "BoxCent= ers" ;ab hier neu einfuegen

res@cnLevelSelectionMode = "Expl= icitLevels"
res@cnLevels = ispan(16,27,1)
res@lbLabelStrings = = (/"ET", "EF", "Dwd", "Dwc", "Dwb", "Dwa", "Dsd", "Dsc", "Dsb", "Dsa", "= Dfd", "Dfc", "Dfb", "Dfa", "Cwc", "Cwb", "Cwa", "Csc", "Csb", "Csa", "Cfc",= "Cfb", "Cfa", "BSh", "BSk", "BWh", "BWk", "Aw", "As", "Am", "Af"/)
= res@cnMinLevelValF = min(Klimaklasse(0,:,:)) - 0.5
res@cnMaxLevelVa= lF = max(Klimaklasse(0,:,:)) + 0.5
res@cnLevelSpacingF = 1.0

res@pmTickMarkDisplayMode = "Always"; use NCL default lat/lon la= bels
res@mpMinLatF      = -23
res@m= pMaxLatF      = -10
res@mpMinLonF &nb= sp;    = -72
res@mpMaxLonF    &n= bsp; = -55
res@mpLimitMode = "latlon"
res@mpLeftCornerLatF= = lat2d(0,0)
res@mpLeftCornerLonF = lon2d(0,0)
res@mpRight= CornerLatF = lat2d(latsize-1,lonsize-1)
res@mpRightCornerLonF = l= on2d(latsize-1,lonsize-1)

res@gsnCenterString = "Klimaklassi= fikation"
plot = gsn_csm_contour_map(wks,Klimaklasse(0,:,:),res)
;res@gsnCenterString = "Mitteltemperatur 1991-2020"
;plo= t(1) = gsn_csm_contour_map(wks,N,res)

;res@gsnCenterString = = "Mitteltemperatur 2021-2050"
;plot(2) = gsn_csm_contour_map(wks= ,O,res)

;res@gsnCenterString = "Mitteltemperatur 2051-2080"<= br /> ;plot(3) = gsn_csm_contour_map(wks,P,res)

;gsn_panel(w= ks,plot,(/2,2/),False)         = ;    ; now draw as one plot

end

I get the following er= ror message:

 

fatal:Subscript out of= range, error in subscript #3
fatal:An error occurred reading Tempera= tur
fatal:An error occurred reading Temperatur
fatal:["Execute.= c":8126]:Execute: Error occurred at or near line 43

 

Line 43 is the followi= ng line: temp = Temperatur(:,0,loc_lat,loc_lon)

 

How can I solve this p= roblem? I hope you can help me.

 

Attachment:

 

- file: MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_Celsiustt_ymonm= ean.nc

- file: MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minus_1_hour_ymonsum_mea= n_30_years.nc

- file:= Koeppen_Geiger_neu.ncl (my NCL code)

- file: ll.out

 

 

Best re= gards

 

Stephan=

_______________________________________________ ncl-talk mailing list =
List instructions, subscriber options, unsubscribe: http://mailman.ucar.edu/mailman/list=
info/ncl-talk 

_______________________________________________
ncl-talk mailin= g list
List instructions, subscriber options, unsubscribe:
http://mailman.uc= ar.edu/mailman/listinfo/ncl-talk