Skip to content

Instantly share code, notes, and snippets.

@pwolfram

pwolfram/KoppenGeiger.m

Created December 13, 2018 16:44
Show Gist options
  • Save pwolfram/376311e23acdc2e1f9159646051ce3f6 to your computer and use it in GitHub Desktop.
Save pwolfram/376311e23acdc2e1f9159646051ce3f6 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
KoppenGeiger.m
function [Class, BroadClass] = KoppenGeiger(T,P) %,classes)
% [Class, BroadClass] = KoppenGeiger(T,P,classes)
%
% Classify monthly temperature and preciptiation climatologies according
% to the Koppen-Geiger climate classification. The T and P inputs
% represent the temperature and preciptiation climatologies,
% respectively. T and P should be provided as three dimensional arrays
% with the third dimension representing time (12 months). The temperature
% data should have units degrees Celsius and the precipitation data units
% mm/month.
%
% The classes input relates the class symbols (e.g., Dwa) to specific and
% broad numeric identifiers (e.g., 21 and 4, respectively). classes
% should be a cell array with the first column containing the numeric
% identifier of the specific class, the second column containing the
% numeric identifier of the broad class, and the third column containing
% the symbol. Other columns are not used in this function. This is an
% example of the classes input:
%
classes = {...
1,1,'Af','Tropical, rainforest',[0 0 255];...
2,1,'Am','Tropical, monsoon',[0 120 255];...
3,1,'Aw','Tropical, savannah',[70 170 250];...
4,2,'BWh','Arid, desert, hot',[255 0 0];...
5,2,'BWk','Arid, desert, cold',[255 150 150];...
6,2,'BSh','Arid, steppe, hot',[245 165 0];...
7,2,'BSk','Arid, steppe, cold',[255 220 100];...
8,3,'Csa','Temperate, dry summer, hot summer',[255 255 0];...
9,3,'Csb','Temperate, dry summer, warm summer',[200 200 0];...
10,3,'Csc','Temperate, dry summer, cold summer',[150 150 0];...
11,3,'Cwa','Temperate, dry winter, hot summer',[150 255 150];...
12,3,'Cwb','Temperate, dry winter, warm summer',[100 200 100];...
13,3,'Cwc','Temperate, dry winter, cold summer',[50 150 50];...
14,3,'Cfa','Temperate, no dry season, hot summer',[200 255 80];...
15,3,'Cfb','Temperate, no dry season, warm summer',[100 255 80];...
16,3,'Cfc','Temperate, no dry season, cold summer',[50 200 0];...
17,4,'Dsa','Cold, dry summer, hot summer',[255 0 255];...
18,4,'Dsb','Cold, dry summer, warm summer',[200 0 200];...
19,4,'Dsc','Cold, dry summer, cold summer',[150 50 150];...
20,4,'Dsd','Cold, dry summer, very cold winter',[150 100 150];...
21,4,'Dwa','Cold, dry winter, hot summer',[170 175 255];...
22,4,'Dwb','Cold, dry winter, warm summer',[90 120 220];...
23,4,'Dwc','Cold, dry winter, cold summer',[75 80 180];...
24,4,'Dwd','Cold, dry winter, very cold winter',[50 0 135];...
25,4,'Dfa','Cold, no dry season, hot summer',[0 255 255];...
26,4,'Dfb','Cold, no dry season, warm summer',[55 200 255];...
27,4,'Dfc','Cold, no dry season, cold summer',[0 125 125];...
28,4,'Dfd','Cold, no dry season, very cold winter',[0 70 95];...
29,5,'ET','Polar, tundra',[178 178 178];...
30,5,'EF','Polar, frost',[102 102 102];...
};
%
% Copyright 2018 Hylke Beck
% www.gloh2o.org
T_ONDJFM = mean(T(:,:,[10 11 12 1 2 3]),3);
T_AMJJAS = mean(T(:,:,[4 5 6 7 8 9]),3);
tmp = T_AMJJAS>T_ONDJFM;
SUM_SEL = false(size(T));
SUM_SEL(:,:,[4 5 6 7 8 9]) = repmat(tmp,1,1,6);
SUM_SEL(:,:,[10 11 12 1 2 3]) = repmat(~tmp,1,1,6);
clear tmp
Pw = sum(P.*single(~SUM_SEL),3);
Ps = sum(P.*single(SUM_SEL),3);
Pdry = min(P,[],3);
%tmp = P; tmp(~SUM_SEL) = NaN; % Too slow
tmp = single(SUM_SEL); tmp(tmp==0) = NaN;
Psdry = nanmin(P.*tmp,[],3);
Pswet = nanmax(P.*tmp,[],3);
%tmp = P; tmp(SUM_SEL) = NaN; % Too slow
tmp = single(~SUM_SEL); tmp(tmp==0) = NaN;
Pwdry = nanmin(P.*tmp,[],3);
Pwwet = nanmax(P.*tmp,[],3);
MAT = mean(T,3);
MAP = sum(P,3);
Tmon10 = sum(T>10,3);
Thot = max(T,[],3);
Tcold = min(T,[],3);
Pthresh = 2*MAT+14;
Pthresh(Pw*2.333>Ps) = 2*MAT(Pw*2.333>Ps);
Pthresh(Ps*2.333>Pw) = 2*MAT(Ps*2.333>Pw)+28;
B = MAP < 10*Pthresh;
BW = B & MAP < 5*Pthresh;
BWh = BW & MAT >= 18;
BWk = BW & MAT < 18;
BS = B & MAP >= 5*Pthresh;
BSh = BS & MAT >= 18;
BSk = BS & MAT < 18;
A = Tcold >= 18 & ~B; % Added "& ~B"
Af = A & Pdry >= 60;
Am = A & ~Af & Pdry >= 100-MAP/25;
Aw = A & ~Af & Pdry < 100-MAP/25;
C = Thot > 10 & Tcold > 0 & Tcold < 18 & ~B;
Cs = C & Psdry<40 & Psdry<Pwwet/3;
Cw = C & Pwdry<Pswet/10;
overlap = Cs & Cw;
Cs(overlap & Ps>Pw) = 0;
Cw(overlap & Ps<=Pw) = 0;
Csa = Cs & Thot >= 22;
Csb = Cs & ~Csa & Tmon10 >= 4;
Csc = Cs & ~Csa & ~Csb & Tmon10>=1 & Tmon10<4;
Cwa = Cw & Thot >= 22;
Cwb = Cw & ~Cwa & Tmon10 >= 4;
Cwc = Cw & ~Cwa & ~Cwb & Tmon10>=1 & Tmon10<4;
Cf = C & ~Cs & ~Cw;
Cfa = Cf & Thot >= 22;
Cfb = Cf & ~Cfa & Tmon10 >= 4;
Cfc = Cf & ~Cfa & ~Cfb & Tmon10>=1 & Tmon10<4;
D = Thot>10 & Tcold<=0 & ~B; % Added "& ~B"
Ds = D & Psdry<40 & Psdry<Pwwet/3;
Dw = D & Pwdry<Pswet/10;
overlap = Ds & Dw;
Ds(overlap & Ps>Pw) = 0;
Dw(overlap & Ps<=Pw) = 0;
Dsa = Ds & Thot>=22;
Dsb = Ds & ~Dsa & Tmon10 >= 4;
Dsd = Ds & ~Dsa & ~Dsb & Tcold<-38;
Dsc = Ds & ~Dsa & ~Dsb & ~Dsd;
Dwa = Dw & Thot>=22;
Dwb = Dw & ~Dwa & Tmon10 >= 4;
Dwd = Dw & ~Dwa & ~Dwb & Tcold<-38;
Dwc = Dw & ~Dwa & ~Dwb & ~Dwd;
Df = D & ~Ds & ~Dw;
Dfa = Df & Thot>=22;
Dfb = Df & ~Dfa & Tmon10 >= 4;
Dfd = Df & ~Dfa & ~Dfb & Tcold<-38;
Dfc = Df & ~Dfa & ~Dfb & ~Dfd;
E = Thot <= 10 & ~B; % Added "& ~B", and replaced "Thot<10" with "Thot<=10"
ET = E & Thot>0;
EF = E & Thot<=0;
Class = zeros(size(T(:,:,1)),'single');
BroadClass = zeros(size(T(:,:,1)),'single');
for cc = 1:length(classes(:,1))
Class(eval(classes{cc,3})) = classes{cc,1};
BroadClass(eval(classes{cc,3})) = classes{cc,2};
end
tmp = single(~isnan(P(:,:,1)+T(:,:,1)));
tmp(tmp==0) = NaN;
Class = Class.*tmp; % Faster than using a mask
BroadClass = BroadClass.*tmp;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment