Friday, November 20, 2015

Matlab: Converison from Arc/Info Grid to GeoTiff Format

Summary

笔记：记录ANUSPLIN LAPGRD输出.grd文件转换为通用文件Geotiff之方法。

Note

The codes in this case is to be oriented to the files associated with Lat/Long coordinate system only. For others, I would code in the future probably. 
 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
% 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

clear;
clc;

% Source Directory
SouDir='./Input';
% Destination Directory
DesDir='./Output';

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

for ii=1:length(files)
    filename=files(ii).name;
    filepath=[SouDir, '/', filename];
    
    % Output file fullpath
    outname=[filename(1:end-3), 'tif'];
    outpath=[DesDir, '/', outname];
    
    % Conversion
    Gonvert2GTiff(filepath, outpath);
    
    disp([num2str(ii), ':', outpath]);
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
34
35
36
37
38
39
40
41
% 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 Gonvert2GTiff(grdpath, gtiffpath);

    % Get the properties of this .grd file
    GeoInfo=GetProperty(grdpath);
    
    % Get the values of this .grd file;
    MatVal=GetValues(grdpath);
    MatVal=reshape(MatVal, GeoInfo.NCOLS, GeoInfo.NROWS);
    MatVal=MatVal';
    
    
    % Spatial extent
    stLong=GeoInfo.XLLCORNER;
    edLong=stLong+GeoInfo.CELLSIZE*GeoInfo.NCOLS;
    stLat=GeoInfo.YLLCORNER;
    edLat=stLat+GeoInfo.CELLSIZE*GeoInfo.NROWS;
    
    % Write a Data Referenced to Geographic Coordinates
    R=georasterref('RasterSize', [GeoInfo.NROWS, GeoInfo.NCOLS], 'Latlim', [stLat, edLat], 'Lonlim', [stLong, edLong], 'ColumnsStartFrom', 'north');
    
    geotiffwrite(gtiffpath, MatVal, R);



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
37
% 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


% 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

References

[1] 实现气象数据的空间分布.
Posted by at
Labels: ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)