Dolphinese: 小功告成

在熬了两夜，整整两天连登陆学校邮箱都顾不上的不断 scripting 和 debug 之后，我总算完成了研究地区的一部分净初级生产力分布图。

以下是完成最后一步所需的 Bourne (Again?) Shell script，在一台连 seq 都没有的 Windows 机器上，MSYS 1.0 和 GRASS 6.3.0CVS 环境下运行通过：


#!/bin/sh
# calculate NPP using CASA model in GRASS

# import GLDAS data and flip it
cd /c/Users/ofnir/images/GLDAS
#rm evapoGLDAS_CLM10_M_WB_A198106.grb.bin
#rm tempeGLDAS_CLM10_M_EB_A198106.grb.bin
ls evapoGLDAS*.bin > filelist2
ls tempeGLDAS*.bin > filelist3
for ((i=1;i<=270;i=i+1)); do
r.in.bin.exe -f input=`head.exe -n $i filelist2 | tail.exe -n 1` output=`head.exe -n $i filelist2 | tail.exe -n 1` bytes=1 north=90N south=60S east=180E west=180W rows=150 cols=360 anull=9.999e+20 --overwrite
g.region.exe rast=`head.exe -n $i filelist2 | tail.exe -n 1`
r.out.ascii.exe input=`head.exe -n $i filelist2 | tail.exe -n 1` output=/c/Users/ofnir/images/GLDAS/test.txt dp=8 null=*
head.exe -n 6 /c/Users/ofnir/images/GLDAS/test.txt > /c/Users/ofnir/images/GLDAS/evapo.txt
for ((j=1;j<=150;j=j+1)); do
tail.exe -n $j /c/Users/ofnir/images/GLDAS/test.txt | head -n 1 >> /c/Users/ofnir/images/GLDAS/evapo.txt
done;
r.in.ascii.exe -f input=/c/Users/ofnir/images/GLDAS/evapo.txt output=`head.exe -n $i filelist2 | tail.exe -n 1` 'mult=1.0 or read from header' 'nv=* or read from header' --overwrite
done;

for ((i=1;i<=270;i=i+1)); do
r.in.bin.exe -f input=`head.exe -n $i filelist3 | tail.exe -n 1` output=`head.exe -n $i filelist3 | tail.exe -n 1` bytes=1 north=90N south=60S east=180E west=180W rows=150 cols=360 anull=9.999e+20 --overwrite
g.region.exe rast=`head.exe -n $i filelist3 | tail.exe -n 1`
r.out.ascii.exe input=`head.exe -n $i filelist3 | tail.exe -n 1` output=/c/Users/ofnir/images/GLDAS/test.txt dp=2 null=*
head.exe -n 6 /c/Users/ofnir/images/GLDAS/test.txt > /c/Users/ofnir/images/GLDAS/tempe.txt
for ((j=1;j<=150;j=j+1)); do
tail.exe -n $j /c/Users/ofnir/images/GLDAS/test.txt | head -n 1 >> /c/Users/ofnir/images/GLDAS/tempe.txt
done;
r.in.ascii.exe -f input=/c/Users/ofnir/images/GLDAS/tempe.txt output=`head.exe -n $i filelist3 | tail.exe -n 1` 'mult=1.0 or read from header' 'nv=* or read from header' --overwrite
done;

for i in 1 3 5 7 9 11; do
# calculate W scalar, part I
echo "Calculating wscalar$i using data of month $(((i+1)/2+6))"
r.mapcalc "wscalar$i"="min(1,0.5+(`head.exe -n $(((i+1)/2)) filelist2 | tail.exe -n 1`*2592000/atpet`printf "%02d" $(((i+1)/2+6))`))"
i=$((i+1))
echo "Calculating wscalar$i using data of month $(((i+1)/2+6))"
r.mapcalc "wscalar$i"="min(1,0.5+(`head.exe -n $(((i+1)/2)) filelist2 | tail.exe -n 1`*2592000/atpet`printf "%02d" $(((i+1)/2+6))`))"
done;

for i in 13 37 61 85 109 133 157 181 205 229 253 277 301 325 349 373 397 421 445 469 493 517; do
# calculate W scalar, part II
for ((j=1;j<=12;j=j+1)); do
echo "Calculating wscalar$i using data of month $j"
r.mapcalc "wscalar$i"="min(1,0.5+(`head.exe -n $(((i+1)/2)) filelist2 | tail.exe -n 1`*2592000/atpet`printf "%02d" $j`))"
i=$((i+1))
echo "Calculating wscalar$i using data of month $j"
r.mapcalc "wscalar$i"="min(1,0.5+(`head.exe -n $(((i+1)/2)) filelist2 | tail.exe -n 1`*2592000/atpet`printf "%02d" $j`))"
i=$((i+1))
done;
done;

for ((i=1;i<=270;i=i+1)); do
# calculate T scalar
amap=`head.exe -n 255 filelist3 | tail.exe -n 1`
bmap=`head.exe -n $i filelist3 | tail.exe -n 1`
echo "Calculating tscalar$((2*i-1)) and tscalar$((2*i))"
r.mapcalc "tscalar$((2*i-1))"="(0.8+0.02*($amap-273.15)-0.05*($amap-273.15)^2)*(1/(1+exp(0.2*(($amap-273.15)-10-($bmap-273.15))))*1/(1+exp(0.3*(-($amap-273.15)-10+($bmap-273.15)))))"
r.mapcalc "tscalar$((2*i))"="(0.8+0.02*($amap-273.15)-0.05*($amap-273.15)^2)*(1/(1+exp(0.2*(($amap-273.15)-10-($bmap-273.15))))*1/(1+exp(0.3*(-($amap-273.15)-10+($bmap-273.15)))))"
done;

# combine everything together
g.region rast=regionFPAR540
for ((i=1;i<=540;i=i+1)); do
amap=regionFPAR$i
bmap=wscalar$i
cmap=tscalar$i
r.mapcalc "NPP$i"="$amap*$bmap*$cmap*`head -n $i /c/Users/ofnir/Desktop/imageprocess/auxillary/tsi.txt | tail -n 1`"
done;

其中前面两个 for 循环是为了把 GLDAS[1. Global Land Data Assimilation Systems The data used in this study were acquired as part of the mission of NASA's Earth Science Division and archived and distributed by the Goddard Earth Sciences (GES) Data and Information Services Center (DISC).] 的二进制格式图像（预先用 wgrib[2. WGRIB is a program to manipulate, inventory and decode GRIB files. ] 从 GRIB 打包格式转换为二进制格式）读入 GRASS。由于 GLDAS CLM10 数据的扫描是从南到北，与 GRASS 的定义相反，我只好用r.out.ascii 将其读出，用很笨的办法重排各行数据，再重新按 ASCII 读入。

用于估算 NPP 的方法是 Carnegie-Ames-Stanford Approach（CASA）模型[3. Field, Christopher B., James T. Randerson, and Carolyn M. Malmstrom. 1995. Global net primary production: Combining ecology and remote sensing. Remote Sensing of Environment 51, no. 1:74-88. ]，其简单形式是$$!NPP(x,t) = APAR(x,t) \cdot \epsilon(x,t)$$ 其中 APAR 是被植物吸收的光合作用活性照射，$$\epsilon(x,t)$$是当时当地植物的光合作用效率。这个公式可以扩展为$$!NPP = S(x,t) \cdot FPAR(x,t) \cdot \epsilon^{*} \cdot T_{1}(x,t)\cdot T_{2}(x,t) \cdot W(x,t)$$ 其中$$S(x,t)$$是表面辐照，定义为总太阳辐照（Total Solar Surface Irradiance）的一半，$$FPAR(x,t)$$ 是光合作用活性照射中可被植物吸收的比例，是 NDVI 的函数，$$T_{1}$$、$$T_{2}$$和$$W$$是与温度和蒸发有关的因子。$$T_{1}$$与 NDVI 峰值时的气温有关，$$T_{2}$$除了与 NDVI 峰值时气温有关外，也与当时当地气温有关。计算公式为$$!T_{1}(x) = 0.8 + 0.02 \times T_{opt}(x) - 0.0005 \times [T_{opt}(x)]^{2}$$ $$!T_{2}(x,t) = C \cdot 1 / {1 + exp[0.2 \cdot (T_{opt}(x) - 10 - T(x,t)]} \cdot 1 / {1 + exp[0.3 \cdot (-T_{opt}(x) - 10 + T(x,t)]}$$。温度信息是从 GLDAS 数据中提取的。这部分计算在代码53--60行实现。

$$W$$ 是在1和$$0.5+\text{Estimated Evapotranspiration}/\text{Potential Evapotranspiration}$$中取最小值，潜在蒸发量是由 UNEP/日本国立环境研究所的 GRID 项目[4. Global Resource Information Database. UNESCO(1987) through GRID.]取得，估计蒸发量取自 GLADS。计算在代码32--51行实现。$$\epsilon*$$是一个可能全局最大值，需要从现场数据估算。我只是做时序比较，因此就不计算它了。

以下是一些计算结果图像，由于数据分辨率和覆盖不同，可以观察到明显的马赛克:D：

Dolphinese

搜索此博客

2008年4月12日星期六

小功告成

没有评论:

发表评论

博客归档