一个快速进行重投影和重采样的小工具(or 外挂)

这篇日志以及这个工具可能会使需要做以下两种工作的同学受益:
(1)需要对大量的高分辨率栅格数据进行投影转换和重采样
(2)使用MOIDS各种产品进行区域、及全球尺度的研究,也就是需要对多个条带的MODIS产品进行拼接、重投影和重采样工作
本程序对于第二种情况有奇效,如果学会使用本程序,将会极大得减少数据预处理时间,使时间成本大大降低。
咳咳,老夫再唠叨一次,如果你是要使用任何MODIS产品进行多年的区域尺度研究,那么这个工具对于你的数据预处理工作将会大有益处,除非你根本就不在意花多少时间。。。

正文:
这个小程序主要适用于以下情境:做大尺度研究,经常需要对大范围、高分辨率栅格数据进行重投影及重采样。通常,在这种情况下,时间成本很高,仅处理单个文件就可能需要较长时间。而我写的这个小程序可以有效解决这一问题,不过,需要指出的是,这个程序只适用于使用“最近邻法”进行重采样的情况,也就是说重采样之后不会改变原始数据的像元值。

程序原理: 在很多情况下,我们使用诸如MODIS,AMSR-E这类遥感产品时,一般都会使用“最近邻法”(Nearest)来进行重投影及重采样。让我们细想一下,其实,这一步的主要时间都用在了计算输出数据的像元位置,而由于使用的插值方法是最邻近法,实际上我们并没有真的插值,我们只是将原始数据按一定投影规则分配到了输出文件的对应栅格里,而理论+实际上,这种“投影规则”对于不同时间的同“批次”数据来说都是一样的,那么,挖掘机就来了,我们为什么每次都要重新计算一次呢?这岂不是在浪费时间!!!所以,这个程序是直接利用 用户提供的“投影规则” 来进行重投影及重采样。说白了,我实际上没有进行任何投影计算,我只是在“搬运”数据,i.e. 我不是在创造真理,我只是真理的搬运工-_-!。所以,我更愿意把这个程序定义为一个“外挂”,而且是加速挂。

实话说,这个程序的原理极其简单,即便是使用IDL、python、R、matlab等脚本语言,也都可以很快实现。我最初就是在课题组里的工作站上使用R语言进行了覆盖整个青藏高原的地温数据的“搬运式”重采样和重投影,效果极佳,因为完全不需要进行常规的重投影和重采样工作,这使得我原本需要一个月才能完成的重投影工作只用了不到1天就全部做完。

但是我后面还要做很多不同数据源的重投影工作,使用R或者其他脚本来做总感觉有点不尽兴,所以我就写了这样一个小工具,考虑到代码执行效率的问题,程序完全使用C编写(编译环境是VS2015RC)。

下面介绍一下这个程序:(我觉得这对我烂到渣的语言表达能力是个巨大挑战)
主界面

1 Continue Reading →

第一个参数化方案(降雪临界温度+融雪基础温度)

前言:本文描述了一个制作极其粗糙的小程序,其主要功能是“基于气象站和MODIS积雪产品估算降雪临界温度和融雪基础温度”,令人遗憾的是,这个垃圾程序是我博士选题的一部分 。 在写这个程序之前还是有点“野心勃勃”的,but,因为害怕影响毕业,所以我并没有选择使用TerraLib这种开源GIS库(那海量的BUG最终把我吓退了)来进行我的工作,而使用Engine则大大降低了实现难度,导致我现在做任何事情都兴味索然。我的博士选题中的一项重要内容是发展一个基于遥感数据的、实用的新参数化方案,目标是估算出在很多水文模型中会被用到的三个物理参数,即:降雪临界温度、融雪基础温度和温度递减率。本文实际上描述的是第一个参数化方案的实施过程。

一、算法:
MODIS积雪产品可以提供每天两景积雪数据,很容易想到,通过逐日比较积雪数据,可以确定“感兴趣区域”在任意相邻的两天是否发生了“降雪”或“融雪”事件。 这里可供选择的判断发生降雪或融雪的指标有:日内/日间积雪变化、日内/日间积雪反照率变化。

二、实施过程:
本程序是在之前写的一个积雪处理程序(http://user.qzone.qq.com/717191727/blog/1382349333)的基础上进一步编写而成的。
如下图
1
Continue Reading →

8种语言性能大比拼:数据处理篇——以nc数据为例

纯属娱乐,新年第一贴,举办一次编程语言的“性能大赛” ,本次比赛的参赛选手有:
Java    C#    C    C++    Fortran    Python    IDL    awk
共8位,We can see,既有像Java、C#这种高富帅选手(标准库很完备,饭来张口,粉丝无数),又包含了C、Fortran这样的纯屌丝选手(靠近底层,事事亲为),还有像Python、IDL这样的科研宠儿,而高大全、伟光正的C++自然不会甘落人后的,最后,awk作为“民科”类选手,也堂而皇之地参赛了,当然了,想也知道这货是成不了黑马的,。言而总之,这场比赛着实令我期待。
Continue Reading →

将并行计算应用于水文模型——以DHSVM为例(二)

终结篇:不务正业了近一个月,终于有点成果了,但是并行计算也就到此为止了,不能再做了,毕业的压力山大。按我本心来说,单纯做水文模拟实在是没劲,因为这就像,大道理谁都懂,可是真正做起来完全不是那么回事儿,那些“道理”在极其复杂的实际情况下未必就是对的,做这些非常不“靠谱”的试验就像是在玩扫雷,即便我最后通关了,要么是我知道这个地方有雷,而有意地避开了,要么就是我“运气好”,没有踩到雷而已。
扯回正题,先亮一下最后的结果:
12
Continue Reading →

将并行计算应用于水文模型——以DHSVM为例(一)

前言: 本文应该是个搞笑帖,但是如果你是做模型的,并且处于菜鸟阶段(这点最重要了, ),迫切希望提高模型的执行速度的话,那么文中提到的方法以及代码也许会对你有些帮助。

背景:说到并行计算,我想对很多人来说早就如雷贯耳了,并行计算目前有两个比较主流的应用手段:MPI , OpenMP (或者 二者配合)。很多人都“歧视”openmp,推崇MPI,认为openmp只不过是几个编译指令而已,很容易掌握,看起来没有MPI复杂,但是我可以说MPI也只不过是一堆消息传递函数而已吗,我个人认为MPI的主要难度在于对原始算法的修改,因为它一般需要一些数理基础。而OpenMP相对更加容易应用,而且实际上一般情况下,只能用OpenMP,这是因为:(1)MPI适用于“分布式多计算机”,也就是由多个计算机连结的工作站集群,我们所(青藏高原研究所)虽然有阳坤老师的一个由10个计算节点组成的工作站,但是由于我只是“下等”用户,每次最多只会分配给我一个节点,申请的节点多于一个就要无限期排队,这个先决条件就迫使我放弃MPI;(2)OpenMP适用于“集中式多处理器”,一台多核笔记本就完全满足条件了,而且,OpenMP的技术规范内容比较简单,易于掌握,可以实现快速应用,这对于我们这些“急于毕业”的人来说尤为重要。

Continue Reading →

利用VIC模型将气象数据日值插为小时值

搭建水文模型,最让人头疼的就是数据不全了,国内能够下载的气象站数据只有温度、相对湿度、风速和降雨,而且还全都是日值。而更为不幸的是,我最近在用的DHSVM模型需要的是小时值,而且还需要短波辐射和长波辐射数据。如何有效利用现有的气象站数据为水文模型做输入就成为我迫切要解决的问题。方法可能有很多,我这里实现了一个解决方案:利用VIC模型来进行插值转换。
题外话:按理说,VIC模型和DHSVM模型与华盛顿大学都渊源很深,可是VIC模型中有相应的插值模块,但是DHSVM中却没有。
VIC的插值算法大体是(可能我的表述不够准确):
(1)首先VIC利用每日平均气温和平均相对湿度估算出每日平均水汽压,在与每日最高温度和最低温度数据以及每日降水数据一起作为输入,使用蒙大拿大学开发的MT-CLIM模型,算出每小时短波辐射 。同时,会算出逐小时云量(后面会用到)
(2)根据上一步算出的每小时短波辐射判断出每日最高温度和最低温度时刻,建立Hermite多项式,保证在最低和最高时刻的一阶导数为0(从而保证了最低温度和最高温度与输入一致),插值得到逐小时的温度。
(3)VIC自身会通过经验公式,利用逐小时的温度和逐日的相对湿度计算逐小时的水汽压,然后利用逐小时的云量、温度和水汽压直接算出逐小时的长波辐射。
以上,是在我详细查看了VIC相关源码的基础上做的一点总结,可能不够准确。
在了解了VIC的工作原理之后,需要做的就是编写相应的接口程序来调用VIC的相应模块生成我们所需要的结果。而且,最好不要修改VIC的源码。所以,我就花了些时间,编写了一个 调用程序(源码贴在后面),没有改动VIC的源文件,只需要将我写的代码文件与VIC的源码文件一起编译即可。
它的使用方法是:
(1)有六个可选项,分别是,”temp”, “prec”, “srad”, “lrad” , “rh”, “wind” ,表示,你希望程序输出温度、降水、短波辐射、长波辐射、相对湿度或者风速。可以多选。如果要计算长波辐射、逐小时温度或者短波辐射,那么需要准备的输入数据是:逐日的降水、相对湿度、最高温度、最低温度和平均温度。如果要计算逐小时的降水,则只需要逐日的降水数据即可。实际上,此程序只会计算短波、长波和温度,其余变量只是简单的平均或填充一天。
(2)此外,用户还需要设置一些必要的参数,程序会在当前目录下寻找一个名为“parameter.txt”的文件,该文件内容如下:

annual_prec    383.96
elevation    4432.4
slope    0
aspect    0
ehoriz    0
whoriz    0
StartYear    2000
StartMonth    1
StartDay    1
latitude    28.58
sw_prec_thresh    0
mtclim_swe_corr     0
have_vp     1
其各项的含义如下,
annual_prec :年平均降水量(mm)
elevation : Site elevation, meters
slope : Site slope, degrees
aspect : Site aspect, degrees (0=N,90=E,180=S,270=W)
ehoriz :Site east horizon, degrees
whoriz :Site west horizon, degrees
StartYear :起始年
StartMonth:起始月
    StartDay:起始日
    latitude :Site latitude, degrees (- for south)
    sw_prec_thresh:Minimum daily precipitation [mm] that can cause “dimming” of incoming shortwave radiation
    mtclim_swe_corr :TRUE = correct MTCLIM’s downward shortwave  radiation estimate in presence of snow
    have_vp :1表示用户提供水汽压数据,0表示使用本程序计算水汽压注意: 通常,slope、aspect、ehoriz以及whoriz这四项是未知的,我查阅了VIC的源码,其在计算时,都是简单地设定为0,如下
1

因此,我在案例中也将这四项全部设为0,至于这四项具体影响有多大。进过简单测试,slope的影响很大,aspect几乎没有影响。slope为0应是假定没有地形影响,以这样的数据再在DHSVM中进行全流域distibute可能会更加合理。

我编写的源码如下:
注意:一个可行的编译运行方法是,下载VIC源码,并解压,将下面的代码复制到名为“vic_invoke.c”的文本文件中,并移动到解压后的目录。Linux下,输入命令:gcc vic_invoke.c Svp.c mtclim_wrapper.c mtclim_vic.c Vicerror.c Close_files.c Nrerror.c Compress_files.c Calc_air_temperature.c Calc_longwave.c -o vic_invoke -Wall -g -I . -lm即可顺利编译,然后./vic_invoke [option]运行。
如果要在windows下运行该程序,可下载安装MinGW,然后将下面源码中的LINUX宏设置为0,然后切换目录到VIC文件夹下,输入同样的编译命令即可生成win32可执行程序。
windows下直接运行,会弹出如下界面
2
Continue Reading →