Monday, November 16, 2015

Matlab: Conversion from Arc/Info Grid to GeoTiff Format Program

Summary

记录ANUSPLIN插值输出结果在投影坐标系下的ASCII文件转换为GeoTiff格式文件。在前面的帖文中，我已经介绍Arc/Info Grid（Geographic Coordinate System）转换为Geotiff代码，此处追加Arc/Info Grid（Projected Coordinate System）转换为Geotiff代码。

Note

代码要求必须有一个Geotiff文件，它包含待转换Arc/Info Grid文件的全部地理信息，这些信息将在代码中赋给转换的矩阵进而输出为Geotiff文件。
Fig. 1
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
% Created by LI Xu
% Version 1.0
% October 26, 2015

% Description
% Main program for conversion to geotiff format 
% from Arc/Info Grid

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

clear;
clc;

% Mask
maskpath='mask.tif';
[~, geo]=geotiffread(maskpath);
info=geotiffinfo(maskpath);

% Source Directory
SouDir='./input';
% Destination Directory
DesDir='./output';

% All files
files=dir([SouDir, '/*.grd']);

for ii=1:length(files)
    filename=files(ii).name;
    % If exist
    otname=strsplit(filename, '.');
    otname=otname{1};
    otpath=[DesDir, '/', otname, '.tif'];
    
    if exist(otpath, 'file')
        disp([num2str(ii), ': ', otname, '.tif']);
        continue;
    end
    
    
    filepath=[SouDir, '/', filename];
    values=GetValues(filepath);
    geoattrs=GetProperty(filepath);
    inmat=reshape(values, geoattrs.NCOLS, geoattrs.NROWS);
    inmat=inmat';
    
    % Export
    
    geotiffwrite(otpath, inmat, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp([num2str(ii), ': ', otname, '.tif']);
end

disp('****************************************************');


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy


function GeoInfo=GetProperty(grdpath)
    

    fid=fopen(grdpath);
    tline=fgets(fid);
    for ii=1:6
        [property, value]=strread(tline, '%s %f');
        str_exe=['GeoInfo.', upper(property{1}), '=', num2str(value), ';'];
        eval(str_exe);
        tline=fgets(fid);
    end
    
    fclose(fid);
    fclose('all');
    
end


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Get values of a .grd file

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

function Output=GetValues(grdpath)
    
    % Read values
    fid=fopen(grdpath);
    
    tline=fgets(fid);
    count=1;
    Output=[];
    while ischar(tline)
        
        if count>6
            Vals=str2num(tline);
            Output=[Output; Vals'];
        end       
        tline=fgets(fid);
        count=count+1;
    end
    
    fclose(fid);
    fclose('all');

end
Attachment.
Posted by at
Labels: ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)