How to use the landmask function in MATLAB?

I want to use the landmask function to exclude the values of the ocean and display the result. However, when I use the landmask function, it says it cannot find 'coast'. Please help.

clc; clear all;

% 파일 이름 설정
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)

% 격자 생성
[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, lat_air, air(:, :, t)', lon_precip, lat_precip, '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(:, :, [1:3, 10:12]), 3); % 겨울 강수량 (144x72)

% 여름/겨울 강수량 비율
summer_ratio = summer_precip ./ MAP; % 여름 강수량 비율
winter_ratio = winter_precip ./ MAP; % 겨울 강수량 비율

% Pthreshold 변수 지정
Pthreshold = zeros(size(MAP)); % 초기화

% 조건 1: 겨울 강수량 비율이 70% 이상인 경우
winter_condition = winter_ratio > 0.7; % 겨울 강수량 조건
Pthreshold(winter_condition) = 2 * MAT(winter_condition);

% 조건 2: 여름 강수량 비율이 70% 이상인 경우
summer_condition = summer_ratio > 0.7; % 여름 강수량 조건
Pthreshold(summer_condition) = 2 * MAT(summer_condition) + 28;

% 조건 3: 나머지 경우
other_condition = ~(winter_condition | summer_condition); % 여름/겨울 외 나머지
Pthreshold(other_condition) = 2 * MAT(other_condition) + 14;

%% Köppen-Geiger 1단계 분류
climate_class = zeros(size(MAP)); % 분류 배열 초기화

% A: Tropical
climate_class(MAP >= 10 & Tcold >= 18) = 1; % 열대 기후

% B: Arid
climate_class(MAP < 10 * Pthreshold) = 2; % 건조 기후

% C: Temperate
climate_class((MAP >= 10) & (Thot > 10) & (Tcold > 0) & (Tcold < 18)) = 3; % 온대 기후

% D: Cold
climate_class((MAP >= 10) & (Thot > 10) & (Tcold <= 0)) = 4; % 냉대 기후

% E: Polar
climate_class((MAP >= 10) & (Thot <= 10)) = 5; % 극지 기후

% 해양 부분 제거
climate_class(land_mask == 0) = NaN; % 해양 부분을 NaN으로 설정

% 지도 시각화
figure;
worldmap('World'); % 전 세계 지도
load coastlines; % 해안선 데이터 로드

% 격자 생성
[X, Y] = meshgrid(lon_precip, lat_precip); % 경도(lon), 위도(lat)를 격자로 변환

% pcolorm 함수 사용
pcolorm(Y, X, climate_class'); % 기후 분류 데이터 시각화
plotm(coastlat, coastlon, 'k', 'LineWidth', 1); % 해안선 그리기

% 지도 설정 추가
c = colorbar; % 컬러바 추가
colormap(parula(5)); % 분류에 맞는 색상 설정
clim([1 5]); % 1~5 범위로 컬러바 설정
c.Ticks = [1, 2, 3, 4, 5]; 
c.TickLabels = {'A','B','C','D', 'E'};

title('Köppen-Geiger Climate Classification (1st Level)'); % 제목 설정

Answer 1

[Angelo Yeo]

As indicated in the landmask's FEX submission, you should fix the following code.

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)'); % 제목 설정

[은진]

Thanks to you, I was able to solve it. I’m so grateful!