Below is the final code along with some fixes.
%%
% https://kr.mathworks.com/matlabcentral/fileexchange/48661-landmask
%%
clc; clear;
warning('off');
% 파일 이름 설정
% fnm = 'air.mon.ltm.1981-2010.nc';
% fnm2 = 'precip.mon.ltm.1981-2010.nc';
%
% % 데이터 읽기
% lon_air = double(ncread(fnm, 'lon')); % 기온 데이터 경도
% lat_air = double(ncread(fnm, 'lat')); % 기온 데이터 위도
% level = ncread(fnm, 'level'); % 고도
% air = ncread(fnm, 'air'); % 기온 데이터 (144x73x12)
%
% lon_precip = double(ncread(fnm2, 'lon')); % 강수량 데이터 경도
% lat_precip = double(ncread(fnm2, 'lat')); % 강수량 데이터 위도
% precip = ncread(fnm2, 'precip'); % 강수량 데이터 (144x72x12)
%
% % 1000 hPa 고도 데이터 선택
% level_idx = find(level == 1000);
% air = squeeze(air(:, :, level_idx, :)); % 1000 hPa 데이터 선택 (144x73x12)
unzip('landmask.zip'); addpath(genpath("./landmask"))
load('dataset.mat')
% 격자 생성
[lon_air_grid, lat_air_grid] = meshgrid(lon_air, lat_air); % 기온 격자
[lon_precip_grid, lat_precip_grid] = meshgrid(lon_precip, lat_precip); % 강수량 격자
% 보간 결과 저장을 위한 배열 초기화
air_interp = zeros(length(lon_precip), length(lat_precip), size(air, 3));
% 각 시간/월별로 보간 적용
for t = 1:size(air, 3)
% air 데이터를 precip의 위경도 격자에 맞춤
air_interp(:, :, t) = interp2(lon_air_grid, lat_air_grid, air(:, :, t)', lon_precip_grid, lat_precip_grid, 'linear')';
end
% air 데이터 보간 완료: air_interp는 precip와 동일한 크기 (144x72x12)
% 연평균 계산
MAT = mean(air_interp, 3); % 연평균 기온 (144x72)
MAP = mean(precip, 3) * 365; % 연평균 강수량 (mm/year) (144x72)
Pdry = min(precip, [], 3); % 가장 건조한 달 강수량 (144x72)
Tcold = min(air_interp, [], 3); % 가장 추운 달의 기온 (144x72)
Thot = max(air_interp, [], 3); % 가장 더운 달의 기온 (144x72)
% 여름 (4월~9월) 강수량 합산
summer_precip = sum(precip(:, :, 4:9), 3); % 여름 강수량 (144x72)
% 겨울 (10월~3월) 강수량 합산
winter_precip = sum(precip(:, :, 10:12), 3); % 겨울 강수량 (144x72)
% 연간 강수량 계산 (144x72)
total_precip = sum(precip, 3); % 전체 강수량 합산 (각 격자에 대해 연간 강수량)
% 여름 강수량 비율 계산
summer_ratio = summer_precip ./total_precip *100; % 여름 강수량 비율 (144x72)
% 겨울 강수량 비율 계산
winter_ratio = winter_precip ./ total_precip*100; % 겨울 강수량 비율 (144x72)
% Pthreshold 변수 초기화
Pthreshold = zeros(size(MAP));
% 여름과 겨울 강수량 비율 계산
% 1. 여름 강수량 비율이 70% 이상인 경우
Psummer = 2 * MAT(summer_ratio > 70) + 28;
% 2. 겨울 강수량 비율이 70% 이상인 경우
Pwinter = 2 * MAT(winter_ratio > 70);
% 3. 나머지 조건 (여름과 겨울 강수량 비율이 모두 70% 미만)
Premainder = 2 * MAT(~(summer_ratio > 70 & winter_ratio > 70)) + 14;
% Pthreshold 계산
for i = 1:size(MAP, 1) % 경도 방향 (144)
for j = 1:size(MAP, 2) % 위도 방향 (72)
if summer_ratio(i, j) > 70
Pthreshold(i, j) = 2 * MAT(i, j) + 28;
elseif winter_ratio(i, j) > 70
Pthreshold(i, j) = 2 * MAT(i, j);
else
Pthreshold(i, j) = 2 * MAT(i, j) + 14;
end
end
end
%% Köppen-Geiger 1단계 분류
climate_class = zeros(size(MAT)); % 분류 배열 초기화
% A: Tropical
climate_class(MAP >= 10*Pthreshold & Tcold >= 18) = 1; % 열대 기후
% B: Arid
climate_class(MAP < 10 * Pthreshold) = 2; % 건조 기후
% C: Temperate
climate_class((MAP >= 10*Pthreshold) & (Thot > 10) & (Tcold > 0) & (Tcold < 18)) = 3; % 온대 기후
% D: Cold
climate_class((MAP >= 10) & (Thot > 10) & (Tcold <= 0)) = 4; % 냉대 기후
% E: Polar
climate_class((MAP >= 10*Pthreshold) & (Thot <= 10)) = 5; % 극지 기후
% 지도 시각화
figure;
worldmap('World'); % 전 세계 지도
load coastlines; % 해안선 데이터 로드
coastlon(coastlon<0)=coastlon(coastlon<0)+360;
% 격자 생성
[lon_precip_grid, lat_precip_grid] = meshgrid(lon_precip, lat_precip); % 강수량 격자 재생성
% 육지 마스크 생성
land_mask = inpolygon(lon_precip_grid, lat_precip_grid, coastlon, coastlat); % 육지 마스크 생성
% figure; imagesc(lat_precip, lon_precip, land_mask)
% geoshow("landareas.shp", "FaceColor","none");
% land_mask 전치 (144x72로 맞추기)
land_mask = double(land_mask');
% 보간: climate_class의 크기와 동일하게 보간
%land_mask_interp = interp2(lon_precip, lat_precip, land_mask, lon_precip, lat_precip, 'linear');
% land_mask_interp의 크기를 climate_class의 크기와 일치시킴
%land_mask_interp = land_mask_interp';
% 육지 마스크가 없는 부분을 NaN으로 설정
% climate_class(land_mask == 0) = NaN;
% pcolorm 함수 사용 시 NaN 처리된 기후 분류 데이터 시각화
h = pcolorm(lat_precip, lon_precip, climate_class'); % 기후 분류 데이터에 육지 마스크 적용
plotm(coastlat, coastlon, 'k', 'LineWidth', 1); % 해안선 그리기
% 컬러맵 설정
load colormap_general4.mat % colrm이 컬러맵 변수라고 가정
% 컬러맵에서 5개의 색상 선택
colrm_5 = colrm(round(linspace(1, size(colrm, 1), 5)), :);
% 컬러맵 뒤집기
colrm_5_reversed = flipud(colrm_5);
% 뒤집힌 컬러맵 적용
colormap(colrm_5_reversed);
% 컬러바 설정
colorbar;
c = colorbar; % 컬러바 추가
clim([1 5]); % 1~5 범위로 컬러바 설정
c.Ticks = [1, 2, 3, 4, 5];
c.TickLabels = {'A','B','C','D', 'E'};
delete(h)
lm = landmask(lat_precip_grid, lon_precip_grid-180); % longitude의 180도 shift
lm = transpose(lm);
lm = circshift(lm, size(lm,1)/2, 1); % longitude의 shift 돌려 놓기
climate_class(~lm) = NaN;
h = pcolorm(lat_precip, lon_precip, climate_class');
title('Köppen-Geiger Climate Classification (1st Level)'); % 제목 설정
