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

Hi Stephan,


take a look at the web page

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

and

  http://www.ncl.ucar.edu/Applications/method_2.shtml

Bye,
Karin

Am 23.06.2014 um 19:57 schrieb "Stephan Herrmann" =
<stephan.w.herrmann@t-online.de>:

> 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 them. 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 =
and 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 =
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_Celsiu=
stt_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_ym=
onsum_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_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 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

_______________________________________________________________________

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


Received on Mon Jun 23 06:12:00 2014