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

From: Karin Meier-Fleischer <meier-fleischer_at_nyahnyahspammersnyahnyah>
Date: Sat Jun 21 2014 - 15:58:16 MDT

Hi Stephan,

the do loops with loc_lat and loc_lon are wrong. You have to change it to

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

To define the variable 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 and plot the climate zones after the K�ppen &
> Geiger climate classification with the following NCL code:
>
> 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"
>
> begin
>
> f =
> addfile("MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_Celsiustt_ymonmean.nc","r")
> Temperatur = f->var167(:,:,:,:)
>
> 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)
>
> 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(:,:,:)
>
> 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)
>
> 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,typeof(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)
> mean_temp = mean_temp + temp(i)
> end do
> mean_temp = mean_temp / 12.0
>
> ;**** Klimaklassifikation 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 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(summe_nied * 0.3 .gt. sommer_nied) then
> regen = 0.0 ; Winterregen
> else
> regen = 7.0 ; gleichm��ige Niederschlagsverteilung
> end if
> end if
>
> warm = 0.0
>
> do k=0,11
> if(temp(k) .gt. 10.0) then
> warm = warm + 1.0
> 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
> Savannenklima ohne Winter
>
> if(min(nied(:)) .gt. 60.0) then
> klasse = klasse + "f" ; Tropisches Regenwaldklima
> else if (min(nied(:)) .gt. 0.04 * (2400.0 - summe_nied)) then
> ;aendern siehe Zettel
> klasse = klasse + "m" ; Mittelform, Regenwaldklima trotz
> einer Trockenzeit
> else
> klasse = klasse + "w" ; 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" ; Steppenklima
> else
> klasse = klasse + "W" ; Wuestenklima
> end if
>
> if(mean_temp .gt. 18.0) then
> klasse = klasse + "h" ; hei�
> else if(mean_temp .le. 18.0) then
> 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" ; Warm-gemae�igtes Klima
>
> 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
>
> else 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
> klasse = klasse + "s" ; warmes sommertrockenes Klima
> else
> klasse = klasse + "f" ; feuchtgemaessigtes Klima
> end if
> end if
>
> 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
> klasse = klasse + "c" ; kuehle Sommer
> end if
> end if
> end if
>
> else if((min(temp(:)) .lt. -3.0) .and. (max(temp(:)) .gt. 10.0)) then
> klasse = "D" ; Boreales oder Schnee-Wald-Klima
>
> 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" ; wintertrockenkaltes Klima
> else
> klasse = klasse + "f" ; winterfeuchtkaltes Klima
> end if
>
> if(max(temp(:)) .gt. 22.0) then
> klasse = klasse + "a" ; hei�e Sommer
> else if((max(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) then
> klasse = klasse + "b" ; warme Sommer
> else if((max(temp(:)) .lt. 22.0) .and. (warm .le. 3.0) .and. (warm
> .ge. 1.0)) then
> klasse = klasse + "c" ; kuehle Sommer
> else if(min(temp(:)) .lt. -38.0) then
> klasse = 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 = klasse + "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
> Klimaklasse(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
> 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
> Klimaklasse(0,loc_lat,loc_lon) = 24.1
> else 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
> Klimaklasse(0,loc_lat,loc_lon) = 21.1
> else if (klasse .eq. "Cfc") then
> Klimaklasse(0,loc_lat,loc_lon) = 20.1
> else if (klasse .eq. "Csa") then
> 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
> 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
> 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,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
> Klimaklasse(0,loc_lat,loc_lon) = 10.1
> else if (klasse .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,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
> 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
> Klimaklasse(0,loc_lat,loc_lon) = 0.1
> 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 if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
>
>
> end do
> end do
>
> wks = gsn_open_wks("ps","Koeppen") ; open a ps file
> ;plot = new(4,graphic) ; create a plot array
>
> gsn_define_colormap(wks,"cb_9step")
>
> res = True
> res@gsnMaximize = True
> res@gsnAddCyclic = False
> res@mpDataBaseVersion = "mediumres"
> res@mpOutlineBoundarySets = "national"
>
> 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 = "BoxCenters" ;ab hier neu einfuegen
>
> res@cnLevelSelectionMode = "ExplicitLevels"
> 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@cnMaxLevelValF = max(Klimaklasse(0,:,:)) + 0.5
> res@cnLevelSpacingF = 1.0
>
> res@pmTickMarkDisplayMode = "Always"; use NCL default lat/lon labels
> res@mpMinLatF = -23
> res@mpMaxLatF = -10
> res@mpMinLonF = -72
> res@mpMaxLonF = -55
> res@mpLimitMode = "latlon"
> res@mpLeftCornerLatF = lat2d(0,0)
> res@mpLeftCornerLonF = lon2d(0,0)
> res@mpRightCornerLatF = lat2d(latsize-1,lonsize-1)
> res@mpRightCornerLonF = lon2d(latsize-1,lonsize-1)
>
> res@gsnCenterString = "Klimaklassifikation"
> plot = gsn_csm_contour_map(wks,Klimaklasse(0,:,:),res)
>
> ;res@gsnCenterString = "Mitteltemperatur 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 = "Mitteltemperatur 2051-2080"
> ;plot(3) = gsn_csm_contour_map(wks,P,res)
>
> ;gsn_panel(wks,plot,(/2,2/),False) ; now draw as one plot
>
> end
>
> I get the following error message:
>
> fatal:Subscript out of range, error in subscript #3
> fatal:An error occurred reading Temperatur
> fatal:An error occurred reading Temperatur
> fatal:["Execute.c":8126]:Execute: Error occurred at or near line 43
>
> Line 43 is the following line: temp = Temperatur(:,0,loc_lat,loc_lon)
>
> How can I solve this problem? I hope you can help me.
>
> Attachment:
>
> - file:
> MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_Celsiustt_ymonmean.nc
>
> - file:
> MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minus_1_hour_ymonsum_mean_30_years.nc
>
> - file: Koeppen_Geiger_neu.ncl (my NCL code)
>
> - file: ll.out
>
> Best regards
>
> Stephan
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk


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


Received on Sat Jun 21 09:58:17 2014

This archive was generated by hypermail 2.1.8 : Wed Jul 23 2014 - 15:33:46 MDT