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

Hi Stephan,

You might take a peek at the instructions at:

    http://ncl.ucar.edu/Applications/o-netcdf.shtml

"Method 1" is pretty straightforward.

Rick



On Mon, Jun 23, 2014 at 11:57 AM, Stephan Herrmann <
stephan.w.herrmann@t-online.de> wrote:

> Hello Dennis, hello Karin,
>
>
>
> thank you very much for your help. It is working now. Now I want to save
> the calculated climate zones in a new NetCDF file instead of plotting the=
m.
> Is that possible and how can I do this?
>
>
>
>
>
> Best regards
>
>
>
> Stephan
>
>
>
>
>
>
>
> -----Original-Nachricht-----
>
> Betreff: Re: [ncl-talk] Calculation and plotting of the climate zones
> after the Köppen & Geiger climate classification
>
> Datum: Mon, 23 Jun 2014 16:09:04 +0200
>
> Von: Dennis Shea <shea@ucar.edu>
>
> An: stephan.w.herrmann@t-online.de
>
>
>
>
> I quickly looked at your code. It looks like you are programming in C
> or fortran? Your are prealloccating memory (malloc, allocate) and the
> looping (for, do). All NCL functions automatically allocate array space a=
nd
> there
> are functions to average and sum. Use of built-in functions makes for
> cleaner and more efficient code. Please read the mini-Language manual at:
> http://www.ncl.ucar.edu/Document/Manuals/
>
> ====
> -->[1]<============
>
> Klimaklasse = new(dimsizes(Niederschlag),float)
> Klimaklasse = 99999.9
>
> can be done in one line
>
> Klimaklasse = new(dimsizes(Niederschlag),float, 99999.9)
>
> -->[2]<============
>
>
> 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
>
> can be replaced by:
>
> For example
>
> mean_temp = dim_avg_n(Temperatur(:,0,loc_lat,loc_lon) , 0)
> summe_nied = dim_sum_n( Niederschlag(:,loc_lat,loc_lon), 0)
>
> or, if you stored the variables in local arrays
>
> mean_temp = dim_avg_n(temp , 0)
> summe_nied = dim_sum_n(nied , 0)
>
> or, if meta data are desired
>
> mean_temp = dim_avg_n_Wrap(Temperatur(:,0,loc_lat,loc_lon), 0)
> summe_nied = dim_sum_n_Wrap(Niederschlag(:,loc_lat,loc_lon), 0)
>
> =======
> ...
> do k=0,11
> if(temp(k) .gt. 10.0) then
> warm = warm + 1.0
> end if
> end do
>
> can be replaced with
>
> warm = dim_num_n(temp.gt.10.0, 0)
>
> ...
>
> 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
>
> replaced by
>
> klasse = where(summe_nied / 10.0 .gt. mean_temp + regen \
> ,klasse + "S", klasse + "W")
> klasse = where(mean_temp .gt. 18.0, klasse + "h", klasse + "k")
>
> etc.
>
> ++++++++
> Please look at the NCL function list.
>
> Good Luck
>
>
>
>
> On Sat, Jun 21, 2014 at 3:58 PM, Karin Meier-Fleischer <
> meier-fleischer@dkrz.de> wrote:
>
>> Hi Stephan,
>>
>> the do loops with loc_lat and loc_lon are wrong. You have to change it t=
o
>>
>> 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 & G=
eiger
>> 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_Cels=
iustt_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, 1=
2),
>> 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 Niederschlagsve=
rteilung
>> 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 Savannenkl=
ima
>> 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 eine=
r
>> 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 hie=
r
>> 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", "A=
w",
>> "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_ymo=
nmean.nc
>>
>> - file:
>> MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minus_1_hour_ymonsum_m=
ean_30_years.nc
>>
>> - file: Koeppen_Geiger_neu.ncl (my NCL code)
>>
>> - file: ll.out
>>
>>
>>
>>
>>
>> Best regards
>>
>>
>>
>> Stephan
>>
>>
>> _______________________________________________ ncl-talk mailing list Li=
st instructions, subscriber options, unsubscribe: http://mailman.ucar.edu/m=
ailman/listinfo/ncl-talk
>>
>>
>> _______________________________________________
>> 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 Mon Jun 23 07:06:02 2014