Matlab制作陆地掩膜

clc;
clear all;
res = 0.05; % grid size of land mask
Latlim = [30 42.5];
Lonlim = [116.5 129.25];
% m_proj('Mercator','lat',[30 42.5],'lon',[116.5 129.25]);
[S,~] = shaperead('landareas');

landmask = false(diff(Latlim)/res,diff(Lonlim)/res);
ref = makerefmat('RasterSize', size(landmask), ...
   'Latlim', Latlim, 'Lonlim',Lonlim ,...
   'ColumnsStartFrom','north','RowsStartFrom','west');

for i=1:length(S)
    n = find(isnan(S(i).X), 1, 'first') - 1;
    [r,c] = map2pix(ref, S(i).X(1:n), S(i).Y(1:n));
    mask = roipoly(landmask,c,r);
    landmask(mask)=true;
end

geotiffwrite('landmask',uint8(landmask),ref)