Wednesday, June 29, 2016

Data: FROM-GLC

Introduction

Global land cover data are key sources of information for understanding the complex interactions between human activities and global change. FROM-GLC (Finer Resolution Observation and Monitoring of Global Land Cover) is the first 30 m resolution global land cover maps produced using Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) data. Our long-term goal in FROM-GLC is to develop a multiple stage approach to mapping global land cover so that the results can better meet the needs of land process modeling and can be easily cross-walked to existing global land cover classification schemes.

Data

下载数据有一点问题,投影信息不正确,坐标信息是正确的。这里可以用MODIS数据投影做一补充。
附件里包含演示数据,下面是一般过程。
打开一个空白的Arcmap文件,里面不包含任何空间信息。依次添加数据h4.tif,再是FROM-GLC_30mTile_h27v05_50km_Majority.tif,再加入Code/ne_10m_admin_0_countries.shp,注意这里顺序不能变化,h4文件必须是第一个加进来的文件,初始化Arcmap空间投影设定,矢量文件用来验证FROM-GLC_30mTile_h27v05_50km_Majority.tif文件位置是否准确。
按照上面的顺序,结果如Fig. 1。
Fig. 1
这个顺序安排的目的在于给FROM-GLC_30mTile_h27v05_50km_Majority.tif补上缺失的投影信息,它本身只剩下坐标信息是正确的。
再讲讲FROM-GLC_30mTile_h27v05_50km_Majority.tif文件的正确输出,右键文件选择Data—>Export Data,注意窗口上选择Data Frame (Current),输出为有正确投影信息及坐标信息的文件。
Fig. 2

References

Saturday, June 18, 2016

Data: Global DEM Data

Introduction

全球范围的DEM数据介绍。

GMTED2010

GMTED2010(Global Multi-resolution Terrain Elevation Data). These data have been collected from a variety of sources using aggregation methods. These datasets are best used for working at the continental scale and with very large regions. 7.5 Arc Sec (225m), 15 Arc Sec (450m), and 30 Arc Sec (1km).数据源主要来自SRTM,占比接近70%。

SRTM

SRTM(Shuttle Radar Topography Mission). 数据在2000年采样。第一版是90m空间分辨率,第二版是30m分辨率。 Downlink: http://viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm

ASTER GDEM2

数据是30m空间分辨率。采样时间似乎是2000年。

AW3D30

AW3D30(ALOS World 3D - 30m). 数据是30m空间分辨率。

ALOS DEM[updated 2020-11-14]

数据是12.5米,一个下载的介绍,这是一个视频

References

Wednesday, June 15, 2016

Matlab+GDAL: ANUSPLIN数据处理套装

Introduction

ANUSPLIN插值过程需要做一些数据格式方面的转换,这里就是为它专门制作一个M代码组合套装,调用GDAL函数。
MAIN_TAG_NODATA_GTIFF.m给原Geotiff文件加上指示Nodata的Tag,方便后续转出为Arc/Grid时有NODATA_value这一配置。
MAIN_SETTING_PRJ_FILE.m带给ANUSPLIN输出的单独.grd文件加上投影信息(当然,他们的投影是一致的),方便随后转出为GeoTiff文件。
MAIN_CONVERT_FORMAT_SOLUTION.m是文件格式转换的解决方案,注意以下转换参数所对应的情况,及时调整,GDAL文件格式列表见参考文献[1]。
1
2
3
tarfot='GTiff';
ext='.tif';
files=dir(fullfile(SouDir, '*.grd'));
File.rar.注意查看.prj文件,这里默认参数是WGS_1984_UTM_Zone_50N.

References

Matlab+GDAL: Add NODATA Tag to Geotiff

Introduction

很多Geotiff文件本身的众多Tag没有包含对于Nodata的设定,这就导致日后的处理过程中不得不针对逐个文件再次定义Nodata,若将Geotiff文件转换为Arc/Grid文件,得到的输出文件是不存在NODATA_value这一设置的。
其实,Geotiff文件加上表示Nodata的Tag也很简单,利用GDAL的gdal_translate命令,演示代码:
 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
% June 15, 2016

% Description:
% Add Nodata Tag to TIFF files

% 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
% http://lixuworld.blogspot.com/

clear;
clc;

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

% All Input Files
files=dir(fullfile(SouDir, '*.tif'));
% Nodata Value
nodata=-99;

% Loop
for ii=1:length(files)
    filename=files(ii).name;
    filepath=fullfile(SouDir, filename);
    otpath=fullfile(DesDir, filename);
    % GDAL Command
    strcmd='gdal_translate -of GTiff -a_nodata ';
    strcmd=[strcmd, num2str(nodata), ' ', filepath, ' ', otpath];
    system(strcmd);
    disp([num2str(ii), ': ', filename]);
end

disp('****************************************');
输出结果前后对比如Fig. 1,经上述代码加上指示Nodata的Tag之后就出现了NODATA_value这一配置,而且新输出文件在ArcMap中打开不需要再次配置就可以直接剔除背景区域,见Fig. 2。
Fig. 1
Fig. 2

References

Tuesday, June 14, 2016

Matlab: Image Data Procedure

[1] 特定目录下,全部影像文件对应像元DN数值:
 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
% Created by LI Xu
% Version 1.0
% June 14, 2016

% Description:
% Recognize pixels

% 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
% http://lixuworld.blogspot.com/


function dns=GetDirDN(TarDir, index)

    dns=[];
    files=dir(fullfile(TarDir, '*.tif'));
    % Loop
    for ii=1:length(files)
        filename=files(ii).name;
        filepath=fullfile(TarDir, filename);
        image=imread(filepath);
        DN=image(index);
        dns=AddData(dns, DN);
    end


end

Tool: 网页的时光机

Introduction

Internet Archive (Wayback Machine) 创建于1996年,总部设在美国的洛杉矶市,至今该网站已经保存超过4870亿份网页(2016年6月14日)。 以学校网站为例介绍它强大的保存库,http://www.bnu.edu.cn/。
(1)打开Internet Archive (Wayback Machine)首页
在页面搜索栏中输入目标网址,如Fig. 1,后单击BROWER HISTORY。
Fig. 1
(2)跳转到缓存历史
Fig. 2
如Fig. 2,红框里指出这一网址在1996年10月19日至2016年6月4日之间共计保存1266次。
(3)查看缓存页面
首先要明确,并不是每一天的网址都一个备份存档。这里有哪天的存档才可以看到哪天的历史页面。还是Fig. 2这张页面,全年每一天的列表在页面下方列出,这年(1996年)对学校网址只有3次保存,十月有两次,十二月有一次,Fig. 3。
Fig. 3
单击十月19日,存档页面见Fig. 4,那个时代的网站还别有一番味道。
Fig. 4
我的科学网博客在这上面也有4次存档,见Fig. 5。
Fig. 5
与Internet Archive (Wayback Machine)相似的网站还有北京大学创办的中国Web信息博物馆,这两个网站在时间存档方面可以互补,后者创建时间相对晚上几年,学校网址的存档最早只到2002年,而我科学网的博客则完全没有存档。

References

Wednesday, June 8, 2016

Matlab: Convert DMS to DD

Introduction

经纬坐标DMS格式转换输出为DD格式。输入文件:
Fig. 1
输出文件:
Fig. 2

Monday, June 6, 2016

IDL: MVC

Introduction

MVC(Maximum-value composite)方法计算得出多个矩阵之间对应位置的最大值,空间数据处理中比较常用,尤其是NDVI。这里介绍两个版本的IDL语言MVC函数,其区别在于计算过程的内存占用情况,若输入N个文件,第一版内存占用为N,第二版则仅有1/N,后者相比前者对内存要求较小,适用于单个文件较大的情况。

第一版

 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
;;Created by LI Xu
;;20 November, 2013
;;Version 1.0

;;Discrption
;;MVC, the maximum-value composite method
;;infullname: string array, consist of N files with full path
;;strSaveDir: string, the dirctory of result
;;strflag: string, flag to month or year , etc..


;+
; :Description:
;    Describe the procedure.
;
; :Params:
;    infullname
;    strSaveDir
;    strflag
;
;
;
; :Author: LiXu
;-
pro MVC, infullname, strSavepath

  
  ;;Read one by one
  ok=Query_Tiff(infullname[0], info)
  Cols=info.dimensions[0]
  Rows=info.dimensions[1]
  N_Dim=size(infullname, /n_elements)
  Pics=DblArr(N_Dim, Cols, Rows)
  for ii=0, N_Dim-1 do Begin
    Pics[ii, *, *]=Read_Tiff(infullname[ii], GEOTIFF=GEOVAR)   
  endfor
  
  ;;Run > 
  Result=Pics[0, *, *]
  ;;Loop >
  for jj=1, N_Dim-1 do Begin
    Result=Result>Pics[jj, *, *]   
  endfor
  
  ;;Write a result tiff
  Write_Tiff, strSavePath, Result, GEOTIFF=GEOVAR, /float
  print, strSavePath+' done!'
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
;; Created by LI Xu
;; Version 1.0
;; June 6, 2016

;; Description:
;; MVC Method

;; 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
;; http://lixuworld.blogspot.com/

pro MVC_TinyMomory, fullnames, stroutpath
  ;; Create a blank Matrix 
  ok=Query_Tiff(fullnames[0], info)
  Cols=info.dimensions[0]
  Rows=info.dimensions[1]
  N_Dim=size(fullnames, /n_elements)
  OutMat=DblArr(Cols, Rows)  
  for ii=0, N_Dim-1 do Begin
    image=Read_Tiff(fullnames[ii], GEOTIFF=GEOVAR)  
    OutMat=OutMat>image  
  endfor
  
  ;;Write a result tiff
  Write_Tiff, stroutpath, OutMat, GEOTIFF=GEOVAR, /float
  print, stroutpath+' done!'


end

Matlab版

 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
% Created by LI Xu
% Version 1.0
% June 21, 2016

% Description:
% MVC Matlab Version

% 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
% http://lixuworld.blogspot.com/

% Note;
% fullpathes: full file pathes in cell
% otpath:   output file path               

function MVC4Image(fullpathes, otpath)
    % First Image
    filepath=fullpathes{1};
    [output, geo]=geotiffread(filepath);
    info=geotiffinfo(filepath);
    [row, col]=size(output);
    
    output=output(:);
    for ii=2:length(fullpathes)
        filepath=fullpathes{ii};
        image=imread(filepath);
        image=image(:);        
        output=[output, image];
        output=max(output, [], 2);   
    end
        
    output=reshape(output, row, col);
    geotiffwrite(otpath, output, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
       
end

References

Thursday, June 2, 2016

Matlab: 区域重叠

Introduction

30 m影像数据重采样至1000 m(以下代称A)可能与其他1000 m影像(以下代称B)在有效区域上不完全匹配,Fig. 1上绿色是共同覆盖区域,差别在于边缘,上方蓝色表示B覆盖而A没有,下方红色表示A覆盖而B没有。
Fig. 1
对于这些边缘的填充处理可以参考:众数插值