Data: History Database of the Global Environment

HYDE presents (gridded) time series of population and land use for the last 12,000 years ! It also presents various other indicators such as GDP, value added, livestock, agricultural areas and yields, private consumption, greenhouse gas emissions and industrial production data, but only for the last century.

The HYDE database is developed under the authority of the Netherlands Environmental Assessment Agency.

Math: The Bivariate Normal Distribution

Introduction

设（X, Y）为二维随机变量，若D（X），D（Y），Cov（X, Y）存在，且D（X），D（Y）都大于0，则有

\[{\rho _{XY}} = \frac{{Cov\left( {X,{\rm{ }}Y} \right)}}{{\sqrt {D\left( X \right)} \cdot \sqrt {D\left( Y \right)} }}\]

ρ_XY为相关系数。

假设随机变量X、Y的相关系数ρ_XY存在。当二者相互独立时，则E（XY）=E（X）E（Y），此时Cov（X, Y）=E（XY）—E（X）E（Y）=0，从而ρ_XY=0，即X、Y不相关。反之，若X、Y不相关，X、Y却不一定相互独立。其实，从“不相关”和“相互独立”的含义来看是明显的，不相关只是就线性关系来说的，而相互独立是就一般关系而言的。X、Y相互独立意味着两个变量之间没有任何关系，而X、Y不相关，仅仅说明X、Y之间无线性关系，但并不排除有非线性关系，比如对数关系、平方关系等等。因此，“不相关”是一个比“相互独立”弱得多的概念。

相关系数

设（X, Y）~N（μ₁, μ₂, σ₁, σ₂, ρ），它的概率密度为

\[\begin{array}{c} f\left( {x,y} \right) = \frac{1}{{2\pi {\sigma _1}{\sigma _2}\sqrt {1 - {\rho ^2}} }}\\ \cdot \exp \left\{ {\frac{{ - 1}}{{2\left( {1 - {\rho ^2}} \right)}}\left[ {\frac{{{{\left( {x - {\mu _1}} \right)}^2}}}{{\sigma _1^2}} - 2\rho \frac{{\left( {x - {\mu _1}} \right)\left( {y - {\mu _2}} \right)}}{{{\sigma _1}{\sigma _2}}} + \frac{{{{\left( {y - {\mu _2}} \right)}^2}}}{{\sigma _2^2}}} \right]} \right\} \end{array}\]

试求X和Y的相关系数。

解：可知

\[E\left( X \right) = {\mu _1},E\left( Y \right) = {\mu _2},D\left( X \right) = \sigma _1^2,D\left( Y \right) = \sigma _2^2\]

而

\[\begin{array}{c} {\rm{Cov}}\left( {X,Y} \right) = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\left( {x - {\mu _1}} \right)\left( {y - {\mu _2}} \right)f\left( {x,y} \right){\rm{d}}} } x{\rm{d}}y\\ \left( {{\rm{Let }}t = \frac{1}{{\sqrt {1 - {\rho ^2}} }}\left( {\frac{{y - {\mu _2}}}{{{\sigma _2}}} - \rho \cdot \frac{{x - {\mu _1}}}{{{\sigma _1}}}} \right),u = \frac{{x - {\mu _1}}}{{{\sigma _1}}}} \right)\\ = \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\left( {{\sigma _1}{\sigma _2}\sqrt {1 - {\rho ^2}} tu + \rho {\sigma _1}{\sigma _2}{u^2}} \right){e^{ - \frac{{{u^2} + {t^2}}}{2}}}{\rm{d}}t{\rm{d}}u} } \\ = \frac{{\rho {\sigma _1}{\sigma _2}}}{{2\pi }}\left( {\int_{ - \infty }^{ + \infty } {{u^2}{e^{ - \frac{{{u^2}}}{2}}}{\rm{d}}u} } \right)\left( {\int_{ - \infty }^{ + \infty } {{e^{ - \frac{{{t^2}}}{2}}}{\rm{d}}t} } \right) + \frac{{{\sigma _1}{\sigma _2}\sqrt {1 - {\rho ^2}} }}{{2\pi }}\left( {\int_{ - \infty }^{ + \infty } {u{e^{ - \frac{{{u^2}}}{2}}}{\rm{d}}u} } \right)\left( {\int_{ - \infty }^{ + \infty } {t{e^{ - \frac{{{t^2}}}{2}}}{\rm{d}}t} } \right)\\ = \frac{{\rho {\sigma _1}{\sigma _2}}}{{2\pi }} \cdot \sqrt {2\pi } \cdot \sqrt {2\pi } \end{array}\]

即有

\[{\rm{Cov}}\left( {X,Y} \right) = \rho {\sigma _1}{\sigma _2}\]

于是

\[{\rho _{XY}} = \frac{{{\rm{Cov}}\left( {X,Y} \right)}}{{\sqrt {D\left( X \right)} \sqrt {D\left( Y \right)} }} = \rho \]

这就是说，二维正态随机变量（X，Y）的概率密度中的参数ρ就是X和Y的相关系数，因而二维正态随机变量的分布完全可由X、Y各自的期望、方差以及他们的相关系数确定。

References

[1] 李正耀, 周德强. 大学数学——概率论与数理统计. 北京: 科学出版社. 2009. Pp: 98-101.

Tips: Defining Contiguity

In spatial autocorrelation analysis some measure of contiguity is required. Contiguity has a rather broad definition depending on the research question, however, most analyses in spatial autocorrelation adhere to a common definition of neighborhood relations. Namely, neighborhood relations are defined as either rooks case, bishops case or queens (kings) case. These are rather simple and intuitive as their names suggest (Fig. 1). Rooks case contiguity is by a neighborhood of 4 locations adjacent to each cell, Bishops only considers the diagonals of the relation and queens or kings case considers a neighborhood of eight cells. These are the most common forms of contiguity used in spatial autocorrelation when considering continuous data in a raster format. Of these three the rooks case is the most commonly used and most programs only will compute this particular case.

However, this type of contiguity is not sufficient for vector data formats or irregularly spaced points. In these cases distances to the four or even n nearest neighbors can be used or distance between a variate X and all neighbors can be used.

Fig. 1

File.rar.

References

[1] Global Spatial Autocorrelation Indices - Moran's I, Geary's C and the General Cross-Product Statistic, from University of Ottawa.

Matlab: 小波分析—时间序列的多时间尺度分析（二）

一则小波分析帖子在我的科学网博客上获得了非常大的关注，至今这一帖子保持着留言询问最多的地位。不过，由于2015年11月的一场浩劫，这个帖子附带的示例代码就都烟消云散了，其后，还是有很多同学留言询问，我在这里把代码再回顾一遍，文字基本不变。

本部分对应——Matlab: 小波分析—时间序列的多时间尺度分析（2）。

Introduction

从原帖第3步计算小波系数开始重新讨论，与原帖不同，这里将尝试调用函数cwt()将最大尺度上限延伸至覆盖整个时间序列，得到小波系数，而后依次绘制、分析多时间尺度下的小波波动情况。

Fig. 1

Fig. 1为小波系数实部等值线图。这里呈现约为40年的主振荡周期，在整个时间序列上出现2个偏多中心和1个偏少中心，分别在1909、1998年和1950年。

Fig. 2

小波系数的模值是不同时间尺度变化周期所对应的能量密度在时间域中分布的反映，系数模值愈大，表明其所对应时段或尺度的周期性就愈强。从图 2可以看出，在降水量演化过程中，118~128年时间尺度模值最大（大于1000），横贯整个时间序列。

Fig. 3

小波系数的模方相当于小波能量谱，它可以分析出不同周期的震荡能量。由Fig. 3可知，118~128年时间尺度的能量最强、周期最显著，占据整个研究时域（1894~2010年）。

Fig. 4

小波方差图中（Fig. 4）存在较为明显的峰值。其中，最大峰值对应着78a的时间尺度，说明78a左右的周期震荡最强，为年降水量变化的第一主周期；55a时间尺度对应着第二峰值，为第二主周期（与原帖相同）；第三、第四峰值分别对应着30a和19a的时间尺度，它们依次为降水量的第三和第四主周期。这说明上述4个周期波动控制着降水量在整个时间域内的变化特征。

值得一提的是，小波方差在100~128a时间尺度上一直处于上升趋势，未出现拐点，这一现象值得再做考虑。

根据小波方差检验的结果，绘制演变的第一和第二主周期小波系数图。

Fig. 5

Fig. 6

从主周期趋势图中可以分析出在不同的时间尺度下，降水量存在的平均周期及丰—枯变化特征。Fig. 5显示，在78a特征时间尺度上，降水量变化的平均周期为50a左右，大约经历了2个丰—枯转换期；而在55a特征时间尺度上（Fig. 6），平均变化周期为35a左右，大约3个周期的丰—枯变化。

附上本帖数据及代码.

More

[1] A Practical Guide to Wavelet Analysis.

[2] 初学应用举例（小波分析，实部画图，小波方差画图及小波模画图等）.

[3] Matlab: 小波分析—时间序列的多时间尺度分析（一）.

Matlab: 小波分析—时间序列的多时间尺度分析（一）

一则小波分析帖子在我的科学网博客上获得了非常大的关注，至今这一帖子保持着留言询问最多的地位。不过，由于2015年11月的一场浩劫，这个帖子附带的示例代码就都烟消云散了，其后，还是有很多同学留言询问，我在这里把代码再回顾一遍，文字基本不变。

本部分对应——MATLAB：小波分析—时间序列的多时间尺度分析。

Introduction

时间序列（Time Series）是地学研究中经常遇到的问题。在时间序列研究中，时域和频域是常用的两种基本形式。其中，时域分析具有时间定位能力，但无法得到关于时间序列变化的更多信息；频域分析（如Fourier变换）虽具有准确的频率定位功能，但仅适合平稳时间序列分析。然而，地学中许多现象（如河川径流、地震波、暴雨、洪水等）随时间的变化往往受到多种因素的综合影响，大都属于非平稳序列，它们不但具有趋势性、周期性等特征，还存在随机性、突变性以及“多时间尺度”结构，具有多层次演变规律。对于这类非平稳时间序列的研究，通常需要某一频段对应的时间信息，或某一时段的频域信息。显然，时域分析和频域分析对此均无能为力。

20世纪80年代初，由Morlet提出的一种具有时—频多分辨功能的小波分析（Wavelet Analysis）为更好的研究时间序列问题提供了可能，它能清晰的揭示出隐藏在时间序列中的多种变化周期，充分反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计。

本帖以美国某气象站1894~2010年连续的年降水量为例，试用小波分析，完成如下任务：①小波变换系数；②绘制小波系数实部等值线图；③绘制小波系数模和模方等值线图；④绘制小波方差图；以及⑤绘制不同时间尺度的小波实部过程线。所谓年降水量时间序列的多时间尺度是指：年降水量在演化过程中，并不存在真正意义上的变化周期，而是其变化周期随着研究尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中。也就是说，年降水量变化在时间域中存在多层次的时间尺度结构和局部变化特征。

小波分析的计算过程请参考：小波分析—经典、小波方差制作步骤等。

小波分析的基本理论在此不多叙述，请参考其他文献。本帖主要介绍小波分析的一般过程，数据及部分代码将附在文末。

0．年降水量的变化趋势分析

该站点的年降水量变化情况（LI_plot函数）如图 0所示，其发展呈微微上升的趋势，降水量最高年份出现在2003年，全年累计达到1610.70mm，最低值出现在1965年，累计仅有748.90mm。

Fig. 0

1．数据的加载

加载多年降水量至MATLAB使用Load命令。

2．边界效应的消除或减小

由于本例中的实测年降水量数据为有限时间数据序列，在时间序列的两端可能会产生“边界效用”。为消除或减小序列开始点和结束点附近的边界效应，须对其两端数据进行延伸。在进行完小波变换后，去掉两端延伸数据的小波变换系数，保留原数据序列时段内的小波系数。本例中，我们利用Matlab小波工具箱中的信号延伸（Signal Extension）功能，对数据两端进行对称性延伸。

具体步骤：在Matlab界面的“Command Window”中输入小波工具箱调用命令“wavemenu”，按Enter键弹“Wavelet Toolbox Main Menu”（小波工具箱主菜单）界面（Fig. 1上）；然后单击“Signal Extension”，打开Signal Extension /Truncation窗口，单击“File”菜单下的“Load Signal”，选择prec.mat文件单击“打开”，出现信号延伸界面。Matlab的Extension Mode菜单下包含多种延伸方式和Direction to extend菜单下的3种延伸模式（Both、Left and Right），在这里我们选择对称性两端延伸进行计算。具体操作过程是：在Extension Mode下选择“Symmetric(Whole-Point)”，Dircetion to extend下选择“Both”，单击“Extend”按钮进行对称性两端延伸计算（Fig. 1下），然后单击“File”菜单下的“Save Tranformed Signal”，将延伸后的数据结果存为ePrec.mat文件。

Fig. 1

从ePrec文件可知（entendvalue函数），系统自动将原时间序列数据向前对称延伸5个单位，向后延伸6个单位。

3．计算小波系数

选择Matlab小波工具箱中的Morlet复小波函数对延伸后的数据序列（ePrec.mat）进行小波变换，计算小波系数并保存。

小波工具箱主菜单界面见Fig. 1上，单击“One-Dimensional”下的子菜单“Complex Continuous Wavelet 1-D”，打开一维复连续小波界面，单击“File”菜单下的“Load Signal”按钮，载入时间序列数据ePrec.mat。Fig. 2第一排的左侧为信号显示区域，右侧区域给出了信号序列和复小波变换的有关信息和参数，主要包括数据长度（Data Size）、小波函数类型（Wavelet：cgau、shan、fbsp和cmor）、取样周期（Sampling Period）、周期设置（Scale Setting）和运行按钮（Analyze），以及显示区域的相关显示设置按钮。本例中，我们选择cmor (1-1.5)、取样周期为1、最大尺度为64（延伸后的时间序列的一半），单击“Analyze”运行按钮，计算小波系数。然后单击“File”菜单下的“Save Coefficients”，保存小波系数为cePrec.mat文件。

Fig. 2

4．计算Morlet复小波系数的实部

去除两端延伸数据的小波系数（entendvalue函数），并计算小波系数实部（real函数）。

5．绘制小波系数实部等值线图

这部分过程应用LI_contourf函数，年降水量小波系数实部等值线图见Fig. 3。

如Fig. 3所示的小波系数实部等值线图。其中，横坐标为时间（年份），纵坐标为时间尺度，图中的等值曲线为小波系数实部值。小波系数实部等值线图能反映年降水量序列不同时间尺度的周期变化及其在时间域中的分布，进而能判断在不同时间尺度上，年降水量的未来变化趋势。

Fig. 3

由Fig. 3可以看出降水量1894~2010年演化过程中存在着28~36a的主振荡周期（这一周期是多个振荡周期叠加的矢量和，随后将在小波方差图中对这一主周期进行分解，剥离出第一主周期等等）。在整个时间尺度上出现4个偏多中心和3个偏少中心，分别为1899、1938、1972、2007年和1919、1956、1990年。

红框处出现4个偏多中心和3个偏少中心，一对临近的偏多和偏少波动就可以被识别为一个周期，所以是28~36a主震荡周期。

周期判断可以从正弦曲线上类比，如下图，两个最邻近波峰（或波谷）之间是一个完整的周期，或者两个最邻近的与x轴向下（或向上）交点也是一个完整的周期。

那么按照正弦曲线的周期检测思路，两个最邻近的偏多中心（或偏少中心）就是一个完整的周期长度。

File.rar.

6．绘制小波系数模和模方等值线图

首先，计算小波系数的模和模方，再利用LI_contourf函数绘制模和模方等值线图。

Fig. 4

Morlet小波系数的模值是不同时间尺度变化周期所对应的能量密度在时间域中分布的反映，系数模值愈大，表明其所对应时段或尺度的周期性就愈强。从Fig. 4可以看出，在降水量演化过程中，纵轴40~64a时间尺度模值最大（大于400），但1924~1994年之间模值小余200，说明在此时段内40~64年时间尺度的周期变化并不明显，1995年之后模值再次增大，此时期后40~64年时间尺度的周期变化趋于显著。

Fig. 5

小波系数的模方相当于小波能量谱，它可以分析出不同周期的震荡能量。由Fig. 5可知，40~64年时间尺度的能量最强、周期最显著，但它的周期变化具有局部性（1924年之前和1994年之后）；10~15年时间尺度能量虽然较弱，但周期分布比较明显，几乎占据整个研究时域（1894~2010年）。

7．绘制小波方差图

小波方差的计算过程参考：小波方差制作步骤（参考文献存在一个错误，小波方差未除N，本贴小波方差计算已有修正）。

小波方差图能反映降水量时间序列的波动能量随尺度（a）的分布情况。可用来确定降水量演化过程中存在的主周期。

Fig. 6

8．主周期趋势图的绘制及其在多时间尺度分析中的作用

根据小波方差检验的结果，绘制演变的第一和第二主周期小波系数图。

Fig. 7

Fig. 8

从主周期趋势图中可以分析出在不同的时间尺度下，降水量存在的平均周期及丰-枯变化特征。Fig. 7显示，在55a特征时间尺度上，降水量变化的平均周期为35a左右，大约经历了3个丰—枯转换期；而在30a特征时间尺度上（Fig. 8），平均变化周期为20a左右，大约6个周期的丰—枯变化。

File.rar.

More

[1] A Practical Guide to Wavelet Analysis.

[2] 初学应用举例（小波分析，实部画图，小波方差画图及小波模画图等）.

Data: The Soil Texture for WRF in China

Introduction

The Weather Research and Forecasting (WRF) Model is a next-generation mesoscale numerical weather prediction system designed for both atmospheric research and operational forecasting needs. The model serves a wide range of meteorological applications across scales from tens of meters to thousands of kilometers. The effort to develop WRF began in the latter part of the 1990's and was a collaborative partnership principally among the National Center for Atmospheric Research (NCAR), the National Oceanic and Atmospheric Administration (represented by the National Centers for Environmental Prediction (NCEP) and the (then) Forecast Systems Laboratory (FSL)), the Air Force Weather Agency (AFWA), the Naval Research Laboratory, the University of Oklahoma, and the Federal Aviation Administration (FAA).

WRF can generate atmospheric simulations using real data (observations, analyses) or idealized conditions. WRF offers operational forecasting a flexible and computationally-efficient platform, while providing recent advances in physics, numerics, and data assimilation contributed by developers across the very broad research community. WRF is currently in operational use at NCEP, AFWA, and other centers.

SRF模型的主要输入数据之一即为土壤质地数据，该数据的含义如表一所示，中文名称选自潘小多等2012。

Table. 1

Matlab操作实现中国WRF土壤质地数据的空间分布，输出设定为GeoTIFF格式，如Fig. 1所示。

Fig. 1

中国WRF土壤质地数据由潘小多老师提供，在此向潘老师致以特别的感谢！

How to Cite This Data

Use is free to all if acknowledgement is made.

潘小多, 李新, 冉有华, 刘超 (2012) 下垫面对 WRF 模式模拟黑河流域区域气候精度影响研究. 高原气象, 31, 657-667.

Pan X, Li X, ShiX, Han X, Luo L, Wang L (2012) Dynamic downscaling of near-surface air temperature at the basin scale using WRF-a case study in the Heihe River Basin, China. Frontiers of Earth Science, 6, 314-323.

WRF.rar.

References

[1] WRF Model Users' Page.

[2] The Weather Research & Forecasting Model.

[3] Weather Research & Forecasting (WRF) Model, ATDD.

[4] 16-category WRF Soil Texture, see Table 3: 16-category Soil Categories.

[5] 潘小多, 李新, 冉有华, 刘超 (2012) 下垫面对 WRF 模式模拟黑河流域区域气候精度影响研究. 高原气象, 31, 657-667.

Matlab: AddData函数

Introduction

数据处理过程中，有时并不能预知记录的个数（行数）。若预定义个数不足或过多，都可能产生不确定的结果。基于以上的需求，自编AddData函数。调用方法示例，如Fig. 1。

Fig. 1

AddData.m

function output=AddData(input, data)

    if isempty(input)
        output(1, :)=data;
    else
        output=input;
        rows=size(input, 1);
        output(rows+1, :)=data;

    end

end

Cell版之AddCell.m

function output=AddCell(input, data)

    if isempty(input)
        output{1}=data;
    else
        output=input;
        num=length(output);
        output{num+1}=data;
    end

end

Char版之AddChars.m

function output=AddChars(input, chararr)

    if isempty(input)
        output=chararr;
    else
        output=input;
        rows=size(output, 1);
        output(rows+1, :)=chararr;
    end
    
end

Band版之AddBand.m，多波段矩阵。

function output=AddBand(input, data)

    if isempty(input)
        output=data;        
    else
        output=input;
        [~, ~, band]=size(input);
        output(:, :, band+1)=data;
    end


end

Tips：OA期刊的风险

自阅读凤凰周刊的文章《OA期刊“中国版面费”争议》之后，我深深的认同部分OA（Open Access）期刊存在着圈钱的可能性。

有鉴于此，国外已经有人关注部分OA期刊的正当性，他是University of Colorado Denver的Jeffrey Beall副教授，他有着20年以上的图书馆学术馆员工作经验，他个人列出了涉嫌圈钱的OA期刊目录可供相关人士参考。

LIST OF STANDALONE JOURNALS 涉嫌期刊

LIST OF PUBLISHERS 涉嫌出版者

Geography：The Fraction of Vegetation Cover

Introduction

植被的动态通常是指在植被组成和结构不发生显著改变情况下的周期性变化，如植被的物候变化、植被覆盖度的变化、初级生产力的变化等。植被的演变通常是指植被组成和结构发生显著改变情况下的非周期性变化或长期变化，主要变现为植被的演替和植被带的移动等。

植被覆盖度的三个半经验公式：

\[f = 1 - {\left( {\frac{{NDV{I_\infty } - NDVI}}{{NDV{I_\infty } - NDV{I_S}}}} \right)^{0.6175}}\]

where NDVI_∞ and NDVI_S are the NDVI for vegetation with infinite LAI and bare soil, respectively.

\[f = {\left( {\frac{{NDVI - NDV{I_S}}}{{NDV{I_\infty } - NDV{I_S}}}} \right)^2}\] \[f = \frac{{NDVI - NDV{I_S}}}{{NDV{I_\infty } - NDV{I_S}}}\]

Reference

[1] 陈效逑, 王恒. 1982~2003年内蒙古植被带和植被覆盖度的时空变化[J]. 地理学报, 2009, 64(1): 84~94.

[2] Zhangyan Jiang, Alfredo R. Huete, Jin Chen, et al.. Analysis of NDVI and scaled difference vegetation index retrievals of vegetation fraction[J]. Remote Sensing of Environment, 2006, 101(3): 366~378.

Data：1 km MODIS-based Maximum Green Vegetation Fraction

These data describe annual maximum green vegetation fraction (MGVF), and are based on 12 years (2001-2012) of Collection 5 MOD13A2 normalized difference vegetation index (NDVI) data. Each map shows MGVF for one year (as well as the average, for all years from 2001-2012), based on the annual maximum NDVI and linear mixing models that describe green vegetation fraction (vs. non vegetated area) for each land cover class for each year. Generation of these data is described in Broxton et al., 2014b. The data has been re-gridded from the MODIS sinusoidal grid to a regular latitude-longitude grid, and the map has 21600x43200 pixels (corresponding to a resolution of 30 arc seconds).

Values range from 0 (corresponding to 0% vegetation cover) to 100 (corresponding to 100% vegetation cover).

Reference

[1] 1 km MODIS-based Maximum Green Vegetation Fraction.

Geography: Slope

How Slope works

For each cell, the Slope tool calculates the maximum rate of change in value from that cell to its neighbors. Basically, the maximum change in elevation over the distance between the cell and its eight neighbors identifies the steepest downhill descent from the cell.

Conceptually, the tool fits a plane to the z-values of a 3 x 3 cell neighborhood around the processing or center cell. The direction the plane faces is the aspect for the processing cell. The lower the slope value, the flatter the terrain; the higher the slope value, the steeper the terrain.

If there is a cell location in the neighborhood with a NoData z-value, the z-value of the center cell will be assigned to the location. At the edge of the raster, at least three cells (outside the raster's extent) will contain NoData as their z-values. These cells will be assigned the center cell's z-value. The result is a flattening of the 3 x 3 plane fitted to these edge cells, which usually leads to a reduction in the slope.

Fig. 1

Fig. 2

The Slope algorithm

Slope is commonly measured in units of degrees, which uses the algorithm:

\[slope\_\deg rees = \frac{{180}}{\pi } \times {\rm{ATAN}}\left( {\sqrt {{{\left( {\frac{{dz}}{{dx}}} \right)}^2} + {{\left( {\frac{{dz}}{{dy}}} \right)}^2}} } \right)\]

The rate of change in the x direction for cell e is calculated with the following algorithm:

\[\frac{{dz}}{{dx}} = \frac{{\left( {c + 2f + i} \right) - \left( {a + 2d + g} \right)}}{{8 \cdot x\_cellsize}}\]

The rate of change in the y direction for cell e is calculated with the following algorithm:

\[\frac{{dz}}{{dy}} = \frac{{\left( {g + 2h + i} \right) - \left( {a + 2b + c} \right)}}{{8 \cdot y\_cellsize}}\]

上面这个算法是计算3 x 3窗口不存在Nodata情况的算法，若存在Nodata或处在边缘位置，算法修正为毗邻像元与中心像元高程差最大的方向的角度θ，若中心像元高程是最大值，则该中心像元坡度设置为0，注意这个修正算法得出的结果可能不正确，所以在计算坡度时应该选择比研究区更大的高程数据，规避边缘效应。

示例数据分别在代码与ArcGIS中计算坡度，明显的差别主要分布在栅格数据的边缘，内部不存在明显的差异，之后裁剪出来所需研究区域即可，已经规避边缘影响。

References

[1] How Slope works, from ArcGIS Resource Center.

Matlab+GDAL: HDF-EOS Translation

(1) HDF-EOS Introduction

HDF is the prescribed format for standard data products that are derived from EOS missions. HDF-EOS (Hierarchical Data Format - Earth Observing System) is a self-describing file format for transfer of various types of data between different machines based upon HDF. HDF-EOS is a standard format to store data collected from EOS satellites: Terra, Aqua and Aura. Two versions of HDF-EOS libraries: HDF-EOS2 based on HDF4 and HDF-EOS5 based on HDF5 are developed.

(2) GeoTIFF Introduction

GeoTIFF is a public domain metadata standard which allows georeferencing information to be embedded within a TIFF file. The potential additional information includes map projection, coordinate systems, ellipsoids, datums, and everything else necessary to establish the exact spatial reference for the file. The GeoTIFF format is fully compliant with TIFF 6.0, so software incapable of reading and interpreting the specialized metadata will still be able to open a GeoTIFF format file.

H4

代码按照关键字（keywds）选择性转换输出影像文件。这里输出NDVI和EVI文件。 H4ToGeoTiff.rar.

SPOT-VEGETATION S10 NDVI (HDF4) 原始数据处理代码。

MODIS SWATH H4文件转出为GeoTiff，注意Swath像元不是正方形，所以转出时需要插值。MODSwathToGeoTiff.rar

H5

出现问题：ERROR 4: Unable to open EPSG support file gcs.csv.

Try setting the GDAL_DATA environment variable to point to the

directory containing EPSG csv files.

Failed to process SRS definition: EPSG:4326，它的解决方法。

官方方案：PROBA-V S1/S10 product Matlab reader.附件。

H5文件因其自定义特征导致各类H5文件难以共用同一解决方案，需要结合数据说明加以判定。这里找到PROBA-V官方开发的Matlab代码，以作演示说明，供参考。 附PROBA-V产品说明文档。

[2019.02.07] 面向PROBA-V的H5文件提取软件，附件。每次按照关键词输出，软件运行需要使用gcs.csv文件，可能首次运行返回错误找不到gcs.csv，需在系统变量新增“GDAL_DATA”，将该文件的路径作为变量值，注意重启CMD和MATLAB。输出各类影像需注意查看数据描述，注意scale_factor（伸缩因子）才能得到最终的结果，如NDVI还有_add_offset更需多加留意。Real=DN*a+b.

References

[1] HDF-EOS Tools and Information Center. This site is dedicated to information about HDF-EOS and about tools available to view or work with HDF-EOS and other NASA HDF files. Programming (C, Fortran, IDL® and MATLAB®) and tool examples to access HDF-EOS and other NASA HDF files are also available. It also hosts an archive of presentations made at an annual HDF/HDF-EOS workshop.

[2] GeoTIFF, from Wikipedia, the free encyclopedia.

[3] HDF Group.

[4] Convert from HDF to GeoTiff.

[5] COMPREHENSIVE EXAMPLES.

[6] Display HDF-EOS5 in ENVI.

[7] ArcGIS Open HDF-EOS5.

[8] How to convert (SWI) HDF5 data to GeoTIFF.

[9] H5 Source: Aura OMI Cloud Data Product-OMCLDO2 (Version 003).

[10] HDF-EOS Translation.

Geography: 月度潜在蒸散量

月度潜在蒸散量（Vicente-Serrano，庄少伟，有删改）：

\[{P_{ei}} = \left\{ {\begin{array}{*{20}{c}} 0&{{T_i} \le 1}\\ {1.6d{{\left( {\frac{{10{T_i}}}{I}} \right)}^a} \times 10}&{1 < {T_i} \le 26.5}\\ {{a_1} + {a_2}{T_i} + {a_3}T_i^2}&{{T_i} > 26.5} \end{array}} \right.\] \[a = 0.49239 + 1.792 \times {10^{ - 2}}I - 7.71 \times {10^{ - 5}}{I^2} + 6.75 \times {10^{ - 7}}{I^3}\] \[I = \sum\limits_{i = 1}^{12} {{I_i}} \] \[{I_i} = \left\{ {\begin{array}{*{20}{c}} 0&{{T_i} \le 0}\\ {{{\left( {\frac{{{T_i}}}{5}} \right)}^{1.514}}}&{{T_i} > 0} \end{array}} \right.\] \[\begin{array}{l} {a_1} = - 415.8547\\ {a_2} = 32.2441\\ {a_3} = - 0.4325 \end{array}\]

式中：d为每月的天数除以30（无量纲）；T_i为第i月平均温度（°C），I_i为第i月加热指数（Heat Index）;I为年加热指数；a、a₁、a₂和a₃均为计算过程中间的系数。

Code.rar.

References

[1] 庄少伟. 基于标准化降水蒸发指数的中国区域干旱化特征分析[D]. 兰州:兰州大学, 2013.

[2] Vicente-Serrano SM,Beguería S, López-Moreno JI. A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index[J]. 2009, 23(7):1696-1718.

Geography: 五期土地利用图

Summary

五期全国土地利用图，覆盖时间有1980年代末期、1995年、2000年、2005年、2010年，这五期土地利用图的分类精度如下：

1980年代末期：95%以上，【刘纪远, 张增祥, 庄大方, 等. 20世纪90年代中国土地利用变化时空特征及其成因分析[J]. 地理研究, 2003, 22(1): 1~12.】

1995年：92%以上，【刘纪远, 张增祥, 庄大方, 等. 20世纪90年代中国土地利用变化时空特征及其成因分析[J]. 地理研究, 2003, 22(1): 1~12.】

2000年：95%以上，【刘纪远, 张增祥, 庄大方, 等. 20世纪90年代中国土地利用变化时空特征及其成因分析[J]. 地理研究, 2003, 22(1): 1~12.】

2005年：95%以上，【刘纪远, 张增祥, 徐新良, 等. 21世纪初中国土地利用变化的空间格局与驱动力分析[J]. 地理学报, 2009, 64(12): 1411~1420.】

2010年：94.3%以上，【刘纪远, 匡文慧, 张增祥, 等. 20世纪80年代末以来中国土地利用变化的基本特征与空间格局. 地理学报, 2014, 69(1): 3~14.】

Computation of the maximum NDVI value

The processing chain of the NDVImax computation contains the following three steps: (1) The NDVI normalized frequency distributions were calculated for each of the vegetation types with equal interval of 0.0001. (2) According to the classification accuracy x, the pixels in the distribution range [(1-x)/2, (1+x)/2] are selected for each of the vegetation types. (3) These selected pixels are then calculated again to a normalized frequency distribution. The 95 percentile of the distribution for tall vegetation types and agriculture is assumed to represent vegetation at full cover and maximum activity with an FPAR value close to 1 (here assumed to be 0.95). The 95% NDVI value of agriculture was used to represent all short vegetation types. The 5% desert value is assumed to represent no vegetation and an FPAR value of 0.001 for all landcover types. PROCESS_NDVI.cdr

土地覆盖分类精度统一为90%，某类型植被NDVI下已经剔除≤0之数值，并排除经MODIS转换而得的NDVI：

A. 第二步和第三步未作向正态分布转换处理的方法5% NDVI 0.04155 95% NDVI 0.61565。File_Untransform.rar.

B. 第二步和第三步过程经过Box-Cox transformation，并设置“'BinWidth', 0.0001”，所得结果为5% NDVI 0.00502 95% NDVI 0.55621。File_boxcoxtransform.rar.

References

[1] ZHU Wenquan, PAN Yaozhong, HE Hao, et al.. Simulation of maximum light use efficiency for some typical vegetation types in China[J]. Chinese Science Bulletin, 2006, 51(4): 457~463.

Geography：蒸散量公式

生物温度（BT）是出现植物营养生长范围内的平均温度，在0~30 °C之间，日均温低于0 °C和高于30 °C者均排除在外。其计算公式为：

\[BT = \frac{{\sum t }}{{365}}{\rm{ or }}BT = \frac{{\sum T }}{{12}}\]

其中BT为年（或月）平均生物温度，t为<30 °C且>0 °C的日均温，T为<30 °C且>0 °C的月均温。

可能蒸散率（Potential Evapotranspiration Rate）是一个综合温度与水分平衡的气候指标，通常用来表征和评价植物的气候控制，可能蒸散的方法有许多种，繁简不一，但都是某种形式的热量（温度与辐射）与水分的函数式，可能蒸散与降水的比率即可能蒸散率，可称为干燥度、干燥指数等。利用Holdridge的方法计算可能蒸散和可能蒸散率。计算公式为：

\[PET = BT \cdot 58.93\] \[PER = \frac{{PET}}{P} = \frac{{BT \cdot 58.93}}{P}\]

P是年降水量。

实际蒸散量之斯来伯公式：

\[E = P\left( {1 - {e^{ - \frac{{PET}}{P}}}} \right)\]

周广胜和张新时以植被表面的二氧化碳通量方程（相当于NPP）与水汽通量方程（相当于蒸发散失）之比确定的植被对水的利用效率为基础，根据地球表面两个公认的平衡方程：水量平衡方程和热量平衡方程从能量与水分对蒸发影响的物理过程出发推导出了联系能量平衡方程和水量平衡方程的区域蒸散模式：

\[E = \frac{{P \cdot {R_n}\left( {{P^2} + {R_n}^2 + P \cdot {R_n}} \right)}}{{\left( {P + {R_n}} \right) \cdot \left( {{P^2} + {R_n}^2} \right)}}\]

式中P为年降水量（mm），R_n为净辐射（mm），且

\[{R_n} = \sqrt {PET \cdot P} \cdot \left( {0.369 + 0.598 \cdot \sqrt {\frac{{PET}}{P}} } \right)\]

File.

References

[1] 李镇清, 刘振国, 陈佐忠, 等. 中国典型草原区气候变化及其对生产力的影响. 草业学报, 2003, 12(1):4~10.

Geography: Net Primary Productivity

File.rar.

Code.