% 国土地理院データを読む
clear
%%【変更必要】%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
place = 'C:/Program Files/landuse2000/tiriin_landuse/LUDATA/';
% データのあるフォルダの名前（CからLUDATAまで全て書く）
meshsize = [100 100]; %NS,WE[10m] 300と400の公約数推奨
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 以下，変更不要
files = ls([place,'/NK*']);
[file_num,n] = size(files);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 【ファイル名の説明】
% ファイル名は、(英数字8文字).(英字3文字)のように構成され、
% 最初の1文字「N」は世界測地系対応であることを示し、
% 2文字目が圏域(首都圏は「S」、中部圏は「C」、近畿圏は「K」)を、
% 3文字目が調査時期を、「_」の後の4桁の数字がファイルコードを
% 表しています。
% 最後の3文字については、「TDU」が土地利用データを、
% 「TDG」が行政区域データであることを表しています。
% 例えば、K5_0022.TDUであれば、
% 近畿圏における1996年（第5時期）のファイルコード0022地区の
% 土地利用データであることがわかります。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 【ファイルのフォーマットの説明】
% 土地利用データファイルは、テキスト形式になっており、
% 平面直角座標系において10m×10mの範囲毎に土地利用分類を
% 2桁の数字にコード化し、3km×4km(南北方向×東西方向)単位を
% 1ファイルとして編集整備されています。
% 10mメッシュデータであるので、1ファイルは南北方向に300等分、
% 東西方向に400等分されており、左(西)から400個の分類コードを
% 並べたものにDOSの改行コードを加えたものを1ラインと呼び、
% このラインを上（北）から順に300ライン並べたものが
% 1ファイルになります。(以下の表を参照して下さい。)
% ┏━━━━━━┯━━━━━┯━━━━┯━━━━┯━━━━┯━━━┯━━━━━┯━━━━━┓
% ┃ファイル　　│ライン番号│データ1 │データ2 │データ3 │　～　│データ399 │データ400 ┃
% ┃コード(4桁) │(3桁) 　　│(2桁) 　│(2桁)　 │(2桁)　 │　　　│(2桁)　 　│(2桁)　 　┃
% ┠──────┼─────┼────┼────┼────┼───┼─────┼─────┨
% ┃(例)0022　　│001　 　　│01　　　│02　　　│03　　　│　～　│16　　　　│17　　　　┃
% ┠──────┼─────┼────┼────┼────┼───┼─────┼─────┨
% ┃　　0022　　│002　　　 │01　　　│02　　　│03　　　│　～　│16　　　　│17　　　　┃
% ┠──────┼─────┼────┼────┼────┼───┼─────┼─────┨
% ┃　　～　　　│　　　　　│　　　　│　　　　│　　　　│　　　│　　　　　│　　　　　┃
% ┠──────┼─────┼────┼────┼────┼───┼─────┼─────┨
% ┃　　0022　　│300　　　 │01　　　│02　　　│03　　　│　～　│16　　　　│17　　　　┃
% ┗━━━━━━┷━━━━━┷━━━━┷━━━━┷━━━━┷━━━┷━━━━━┷━━━━━┛
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ns = 300;
we = 400;
x  = fix(300/meshsize(1));
y  = fix(400/meshsize(2));
tic 
output = zeros(x*34,y*25);
for i = 1:file_num
    fid = fopen([place,'/',files(i,:)]);
    for j = 1:300
        head = fscanf(fid,'%4d%3d',[2 1]);
        numdata(j,:) = fscanf(fid,'%2d',[400 1]);
    end
    fclose(fid);
    for j = 1:x
        for k = 1:y
            for l = 1:18
                a = find(numdata(meshsize(1)*(j-1)+1:meshsize(1)*j,...
                         meshsize(2)*(k-1)+1:meshsize(2)*k) == l);
                a(1) = 0;
                b(l) = a(end);
                clear a
            end
            clas(j,k) = find(b == max(b));
        end
    end
    pla = [str2double(files(i,5:6))+1,str2double(files(i,7:8))-4];
    output(x*(pla(1)-1)+1:x*pla(1),y*(pla(2)-1)+1:y*pla(2)) = clas;
    i/398
end
toc
save('landuse','output')
a = 1;
for i = 1:x*34
    for j = 1:y*25
        a = find(output(i,j) == 0);
        if a == 1
            output(i,j) = NaN;
        end
    end
end
pcolor(output)
axis ij equal tight
title([num2str(meshsize(1)/100),'km x ',num2str(meshsize(2)/100),'km'])
