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/nclscri=
pts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contri=
buted.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
f = addfile("MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1=
961-1990_167_minus_1_hour_Celsiustt_ymonmean.nc","r")
Temperatur = f=
->var167(:,:,:,:)
latsize = getfilevardimsizes(f,"rlat")
rotpole = f->rotated_=
pole ; need attributes for correc=
t remap
ll = asciiread("ll.out",(/latsize,lonsize,2/),"float")=
lat2d = ll(::,:,0)
lon2d = ll(::,:,1)
Temperatur=
@lat2d = lat2d
Temperatur@lon2d = lon2d
f = addfile("=
MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minus_1_hour_ymonsum_mean=
_30_years.nc","r")
Niederschlag = f->var142(:,:,:)
lat=
size = getfilevardimsizes(f,"rlat")
lonsize = getfilevardimsizes(f=
,"rlon")
rotpole = f->rotated_pole &=
nbsp; ; need attributes for correct remap
ll = asciiread=
("ll.out",(/latsize,lonsize,2/),"float")
lat2d = ll(::,:,0)
lo=
n2d = ll(::,:,1)
Niederschlag@lat2d = lat2d
Niederschla=
g@lon2d = lon2d
Klimaklasse = new(dimsizes(Niederschlag),flo=
at)
Klimaklasse = 99999.9
do loc=
_lat = 0,max(dimsizes(lat2d))-1
do loc_lon = 0,max(dimsizes(lon2d)=
)-1
temp = new(12,typeof(Temperatur))
temp=
= Temperatur(:,0,loc_lat,loc_lon)
nied = new(12,type=
of(Niederschlag))
nied = Niederschlag(:,loc_lat,loc_lon)
summe_nied = 0.0
mean_temp = 0.0
=
do i = 0,11
summe_nied = summe_nied + nied(i)
end do=
mean_temp = mean_temp / 12.0
;**** Klimaklassifik=
ation nach Koeppen-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 b=
is Oktober?
sommer_nied = sommer_nied + nied(k)
if(summe_nied * 0.7 .lt. sommer_nied) the=
n
regen = 14.0 =
; ; Sommerregen
else i=
f(summe_nied * 0.3 .gt. sommer_nied) then
regen ==
0.0 =
; ; Winterregen
else
regen =
= 7.0 &n=
bsp; ; gleichmäßige Niederschlagsverteilung
&nb=
sp;end if
end if
warm = 0.0
do k=0,11
if(temp(k) .gt. 10.0) then
=
warm = warm + 1.0
&=
nbsp; 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- oder Savannenk=
lima ohne Winter
if(min(nied(:)) .gt. 60.0) t=
hen
klasse = klasse + "f" =
; Tropisches Regenwaldklima
else if (min(nie=
d(:)) .gt. 0.04 * (2400.0 - summe_nied)) then ;aendern siehe Zettel
&n=
bsp; klasse = klasse + "m" ; Mi=
ttelform, Regenwaldklima trotz einer Trockenzeit
el=
se
klasse = klasse + "w" &=
nbsp; ; Savannenklima
end if
=
; end if
else if(summe_nied / 10.0 .lt. 2.0 * (mean_temp +=
regen)) then
klasse = "B" =
; ; Trockenklima
if(summe_nied / 10.0 .gt. mean_temp + regen) then
=
klasse = klasse + "S" &nb=
sp; ; Steppenklima
else
&nbs=
p; klasse = klasse + "W" ; Wuestenklima
&nbs=
p; end if
if(m=
ean_temp .gt. 18.0) then
klasse =
= klasse + "h" ; heiß
else if(m=
ean_temp .le. 18.0) then
klasse =
= klasse + "k" ; kalt
end if
&nb=
sp; end if
else if((min(temp(:)) .gt. -3.0) .and. (min(te=
mp(:)) .lt. 18.0) .and. (summe_nied / 10.0 .gt. 2.0 * (mean_temp + regen)))=
then
klasse = "C" &=
nbsp; ; Warm-gemaeßig=
tes Klima
if(min((/nied(mod(11 + hk, 12)),nie=
d(mod(0 + hk, 12)),nied(mod(1 + hk, 12))/)) .lt. 0.1 * max((/nied(mod(5 + h=
k, 12)),nied(mod(6 + hk, 12)),nied(mod(7 + hk, 12))/))) then
&nb=
sp;
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
&nb=
sp; klasse = klasse + "s" ; warmes sommertrockenes Klim=
a
else
kl=
asse = klasse + "f" ; feuchtgemaessigtes Klima
&nb=
sp; end if
end if
&nbs=
p; if(max(temp(:)) .gt. 22.0) then
=
; klasse = klasse + "a" ;heisse Sommer
=
; else if((max(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) then ;ab hier neu=
einfuegen
klasse = klasse + "b=
" ; warme Sommer
else if((max(temp(:)) =
.lt. 22.0) .and. (warm .le. 3.0) .and. (warm .ge. 1.0)) then
&nb=
sp; klasse = klasse + "c" ; kuehle So=
mmer
end if
end if
&nb=
sp; end if
else if((min(temp(:)) .lt. -3.0) .=
and. (max(temp(:)) .gt. 10.0)) then
klasse = "D"&=
nbsp; &nbs=
p; ; Boreales oder Schnee-Wald-Klima
&n=
bsp; 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
&n=
bsp; klasse = klasse + "w" ; wintertrockenkaltes Klima
&=
nbsp; else
klasse ==
klasse + "f" ; winterfeuchtkaltes Klima
=
; end if
if(max(temp(:)) .=
gt. 22.0) then
klasse = klasse =
+ "a" ; heiße Sommer
else if((ma=
x(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) then
&n=
bsp; klasse = klasse + "b" ; warme Sommer
&n=
bsp; else if((max(temp(:)) .lt. 22.0) .and. (warm .le. 3.0) .an=
d. (warm .ge. 1.0)) then
klasse =
= klasse + "c" ; kuehle Sommer
else i=
f(min(temp(:)) .lt. -38.0) then
k=
lasse = klasse + "d" ; strenge Winter
=
end if
end if
end if
=
end if
else if(max(=
temp(:)) .lt. 10.0) then
klasse = "E" =
; =
Schneeklima
if(max(temp(:)=
) .lt. 0.0) then
klasse = klass=
e + "F" ; Tundrenklima
else
=
klasse = klasse + "T" ; Klima =
ewigen Frostes
end if
end if
end if
end if
end if
end if
if (klasse .eq. "Af") then
&nb=
sp; Klimaklasse(0,loc_lat,loc_lon) = 30.1
else if=
(klasse .eq. "Am") then
Klimaklasse(0,loc_la=
t,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 (klasse .eq. "BWh") then
Klimaklasse(=
0,loc_lat,loc_lon) = 25.1
else if (klasse .eq. "BSk") then
&=
nbsp; else if (klasse .eq. "BSh") then
Klimak=
lasse(0,loc_lat,loc_lon) = 23.1
else if (klasse .eq. "Cfa") t=
hen
Klimaklasse(0,loc_lat,loc_lon) = 22.1
K=
limaklasse(0,loc_lat,loc_lon) = 21.1
else if (klasse .eq. "Cf=
c") then
Klimaklasse(0,loc_lat,loc_lon) = 2=
0.1
else if (klasse .eq. "Csa") then
&nb=
sp; Klimaklasse(0,loc_lat,loc_lon) = 19.1
else if (klasse .eq=
. "Csb") then
Klimaklasse(0,loc_lat,loc_lon) =
= 18.1
else if (klasse .eq. "Csc") then
&nbs=
p; Klimaklasse(0,loc_lat,loc_lon) = 17.1
else if (klass=
e .eq. "Cwa") then
Klimaklasse(0,loc_lat,loc_=
lon) = 16.1
else if (klasse .eq. "Cwb") then
=
; Klimaklasse(0,loc_lat,loc_lon) = 15.1
else if (=
klasse .eq. "Cwc") then
Klimaklasse(0,loc_lat=
,loc_lon) = 14.1
else if (klasse .eq. "Dfa") then
=
Klimaklasse(0,loc_lat,loc_lon) = 13.1
else=
if (klasse .eq. "Dfb") then
Klimaklasse(0,lo=
c_lat,loc_lon) = 12.1
else if (klasse .eq. "Dfc") then
&=
nbsp; Klimaklasse(0,loc_lat,loc_lon) = 11.1
=
else if (klasse .eq. "Dfd") then
Klimaklasse=
(0,loc_lat,loc_lon) = 10.1
else if (klasse .eq. "Dsa") then
&n=
bsp; else if (klasse .eq. "Dsb") then
Klimakl=
asse(0,loc_lat,loc_lon) = 8.1
else if (klasse .eq. "Dsc") the=
n
Klimaklasse(0,loc_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
Klimaklasse(0,loc_lat,loc_lon) ==
5.1
else if (klasse .eq. "Dwb") then
&n=
bsp; Klimaklasse(0,loc_lat,loc_lon) = 4.1
else if (klasse .eq=
. "Dwc") then
Klimaklasse(0,loc_lat,loc_lon) =
= 3.1
else if (klasse .eq. "Dwd") then
=
; Klimaklasse(0,loc_lat,loc_lon) = 2.1
else if (klasse =
.eq. "EF") then
Klimaklasse(0,loc_lat,loc_lon=
) = 1.1
else if (klasse .eq. "ET") then
&nbs=
p; Klimaklasse(0,loc_lat,loc_lon) = 0.1
end if
&nb=
sp; end if
end if
end if
end if
&n=
bsp; end if
end if
end if
end if
&=
nbsp; end if
end if
end if
end if
=
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end if
end do
en=
d do
wks = gsn_open_wks("ps","Koeppen") =
; ; open a ps f=
ile
;plot = new(4,graphic) =
&nb=
sp; ; create a plot array
gsn_defi=
ne_colormap(wks,"cb_9step")
res = True
res@gsnMaximize =
= True
res@gsnAddCyclic = False
res@mpDataBaseVersion = "me=
diumres"
res@mpOutlineBoundarySets = "national"
&nb=
sp;
res@mpCenterLonF = rotpole@=
grid_north_pole_longitude
res@mpCenterLatF = rotpole@gri=
d_north_pole_latitude -90
res@cnFillOn = True
res@cnLinesOn ==
False
res@cnLineLabelsOn = False
res@cnExplicitLegendLab=
elsOn = True
;res@lbLabelAlignment = "BoxCenters" ;ab hier neu ein=
fuegen
res@cnLevelSelectionMode = "ExplicitLevels"
res@cn=
Levels = 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", "BS=
k", "BWh", "BWk", "Aw", "As", "Am", "Af"/)
res@cnMinLevelValF = min=
(Klimaklasse(0,:,:)) - 0.5
res@cnMaxLevelValF = max(Klimaklasse(0,:,=
:)) + 0.5
res@cnLevelSpacingF = 1.0
res@pmTickMarkDisplay=
Mode = "Always"; use NCL default lat/lon labels
res@mpMinLatF &=
nbsp; = -23
res@mpMaxLatF =
= -10
res@mpMinLonF = -72
res@mpMaxLonF = -55
res@mpLimitMode=
= "latlon"
res@mpLeftCornerLatF = lat2d(0,0)
res@mpLeftCorne=
rLonF = lon2d(0,0)
res@mpRightCornerLatF = lat2d(latsize-1,lonsize=
-1)
res@mpRightCornerLonF = lon2d(latsize-1,lonsize-1)
re=
s@gsnCenterString = "Klimaklassifikation"
plot = gsn_csm_contour_m=
ap(wks,Klimaklasse(0,:,:),res)
;res@gsnCenterString = "Mittelt=
emperatur 1991-2020"
;plot(1) = gsn_csm_contour_map(wks,N,res)
=
;res@gsnCenterString = "Mitteltemperatur 2021-2050"
;plot(2) =
= gsn_csm_contour_map(wks,O,res)
;res@gsnCenterString = "Mit=
teltemperatur 2051-2080"
;plot(3) = gsn_csm_contour_map(wks,P,res)
;gsn_panel(wks,plot,(/2,2/),False) &=
nbsp; ; 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 Temperat=
ur
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)
-
Best re= gards
Stephan=
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Wed Jul 23 2014 - 15:33:46 MDT