基于GDAL库的温度推算模型的研究与实现

刘文毫，卫建国,2，张春梅，杨 豫，刘兆宇

(1.北方民族大学 计算机科学与工程学院，宁夏 银川 750021；2.中国气象局旱区特色农业气象灾害监测预警与风险管理重点实验室，宁夏 银川 750002；3.宁夏大学 农学院，宁夏 银川 750021)

0 引 言

为对中宁太阳梁试验区农作物霜冻情况进行研究，需要对研究区域的温度分布状况进行推算。目前对特定区域的温度分布状况进行推算研究的成果较多，Doug M. Smith等[1]通过建立的模型对全球的温度进行了推算，Y.Radhika等[2]采用支持向量机技术对区域的大气温度进行了预测。国内学者在温度推算方面也有很深入的研究。例如，唐圣钧等[3]对贵州省南北两个区域，结合多元回归等方法，分区建立了基于纬度、海拔和坡向以及坡度因子的气温降水模型，在研究区域进行气候要素的小网格推算，结果表明纬度因子对温度和降水的贡献最大。卫建国[4]根据温度分布和热量理论，在250 m的地理网格资料中，采用多点推算方法，综合改进距离加权平均法，对运算中的多层推算结果进行权值订正，结果表明经过多点推算的温度图像平滑而富有变化，从根本上消除了单点推算的台阶现象。刘静等[5]利用线性内插确定界限温度的初终日期，利用逐旬气候资料进行小网格订正，从而得出了红寺堡等地区的温度分布。但是受限于自然条件及其他方面的限制，这些对温度推算的研究都是采用传统实测数据对地区的温度进行分布，而没有采用地形数据结合经纬度信息等资料对区域的温度进行推算，在实现方式上大多是使用传统的C/C++编程技术进行实现，没有采用最新的Python开发技术进行实现，因此处理效率相对较低，无法应对基于大批量数据的温度推算工作。基于此现状，文中基于GDAL开源库，以无人机航测5 m分辨率DEM数据及实测的经纬度数据为数据源，实现了基于小区域温度推算模型的温度分布推算。

1 试验区概况

太阳梁位于卫宁平原之上的中宁县渠口农场，南与中宁县枣园乡接壤，以白马梁为界；北与青铜峡与中宁交界线的胡子沟为界；西倚内蒙古阿拉善左旗头道湖，东靠包兰铁路。试验区地形复杂多变，由土石山地和缓坡丘陵组成，绝大部分为缓坡丘陵，局地为微坡地。

试验区处于中温带干旱区，大陆性季风气候特征明显，光照充足、热量丰富、蒸发强烈、昼夜温差大，气候干旱少雨，春暖快、夏热短、秋凉早、冬寒长。主要气象灾害有干旱、沙尘暴、霜冻、冰雹、暴雨、干热风等。

由于研究区域气候特征独特，农业灾害类型多样，因此该试验区适合农业气象灾害的研究。又因试验区没有详细的温度分布数据，因此，需要建立一套合理的温度推算模型，推算该试验区温度分布状况。

太阳梁桃园试验区布局如图1所示。

2 数据来源及温度推算模型

项目采用太阳梁试验区无人机航测的5 m分辨率细网格DEM数据作为数据源。通常情况下，某地区的平均温度分布情况与该地区大气环境、海拔高度、经纬度及小气候环境有关[6]。不同季节太阳辐射[7]和大气环流背景是导致气温季节变化的主要原因，同时，坡度坡向也会影响该地区的辐射分布情况[8]。该项目采用前人[3]对贵州全区气候要素的推算模型为基础，并结合太阳梁试验区的具体气候特征，设定每项气候要素的权重系数，使其能够更好地实现对太阳梁试验区的温度推算。具体推算模型如下：

图1 太阳梁桃园试验区布局 (左侧及底部阴影区为研究区域)

设基准站的平均气温为S，则预测点的温度值可表示为：

S=f(λ,φ,h,χ,ω)+T

(1)

其中，S为某预测点气候要素值；λ为预测点与基准点的经度差值；φ为预测点与基准点的纬度差值；h为预测点与基准点的海拔差值；χ为预测点与基准点的坡向差值；ω为预测点与基准点的坡度差值。

基于气象科研人员对模型中的每个因子进行分析评估得到各个因子的系数，同时结合研究区域2018年3月的基准点实测的平均温度为14.8℃，确定最终的温度推算模型如下：

T=[1λ+1φ+0.03h+(-0.01)χ+0.01ω]*

0.05+14.8

(2)

3 关键技术

3.1 GDAL地理空间数据抽象库

GDAL[9]全称为地理空间数据抽象库(geospatial data abstraction library)，是众多操作地理空间栅格数据的开源库之一。它利用抽象数据模型来表达所支持的各种文件格式，并使用一系列命令行工具来进行数据的转换和处理。GDAL包含两个库：用于操作地理空间栅格数据[10]的GDAL和用于操作地理空间矢量数据的OGR，且其支持的数据格式也相当广泛，包括GeoTiff(.tiff)、Erdas Image(.img)、(.grd)等155种栅格数据，以及(.shp)、(.gml)、(.GeoPackage)等95种矢量数据。GDAL库还提供了一系列算法的接口[11]，如矢量栅格化、栅格矢量化、遥感数据的空间纠正等，并对这些算法提供了可以运行的文件，方便用户的使用。利用GDAL和OGR这两个开源库，可实现对栅格和矢量数据的操作及处理。该项目采用开源GDAL库读取了.grd格式的原始DEM文件并转换成通用GeoTiff文件，基于GeoTiff文件提取温度推算所需坡度、坡向等因子并写入到栅格文件中。通过调用GDAL库操作栅格数据的具体流程如图2所示。

图2 GDAL库操作栅格数据流程

3.2 Numpy库及Matplotlib库

Numpy库是Python进行科学计算的基础模块，它是一个提供多维数据对象的Python库，包含了多种衍生对象以及一系列为快速计算数组而生的例程[12]。Numpy库最核心的部分是ndarray对象，通过ndarray对象可以实现对栅格数据的有效存储。同时，Numpy库具有强大的矩阵运算能力，可以实现矩阵求逆、求特征值、求特征向量等操作。Numpy科学运算库为温度推算工作的实现提供了支撑。

Matplotlib库[13]是Python编程语言及Numpy扩展包的可视化操作界面，其提供了面向对象的API用于通用GUI工具包。Matplotlib库类似于MATLAB中的Pyplot模块。基于Matplotlib库可将温度推算中的每个因子以图像的形式直观显示，并自动生成图例，生成的图例不仅方便检查数据的准确性，而且可以帮助灵活地调整权重系数，从而实现对温度分布状况更加准确的推算。

4 项目实现流程

项目首先采用无人机航拍的细网格DEM数据生成海拔格点数据，生成该地区的坡度格点数据、坡向格点数据。依据实测的研究区域经纬度范围，生成包含该区域每个格点的经度纬度值的格点数据。由于生成的格点数据结果为标准矩形矩阵，因此需要使用DEM数据生成二值矩阵将格点数据裁剪为研究区形状。获取了温度推算所需5个因子的格点数据以后，基于Python的ConfigParser模块，读取了每个因子的权重系数并基于Numpy库实现模型的运算操作。温度推算技术路线图如图3所示。

图3 温度推算技术路线图

4.1 温度推算因子的提取

(1)海拔因子。

由于DEM数据[14-15]已经包含研究区域的海拔数据，因此只需对DEM数据进行格式转换，将原始.grd格式文件转换为模型推算所需GeoTiff格式。具体步骤如下：通过GDAL.open()函数打开.grd格式的DEM文件，读取为Numpy库支持的ndarray格式，然后通过GDAL库提供的WriteArray()函数将数据写入为GeoTiff格式。

(2)坡度坡向因子。

由于GDAL库包含生成坡度和坡向的函数，因此该项目通过调用GDAL库的函数生成了包含坡度坡向两个因子的文件。坡度坡向因子生成步骤如图4所示。

图4 坡度坡向因子生成步骤

生成坡度坡向因子需调用subprocess模块，subprocess模块可以通过命令行的方式执行GDAL命令，这是执行生成坡度坡向因子命令的前提。引用完成subprocess后开始编写生成坡度坡向因子的命令，需调用GDAL库中gdaldem模块，然后指定如下参数：指定需要生成的因子(slope、aspect)；指定输出图像的压缩方式，该项目采用LZW压缩方式；指定其他参数，如-s指定计算比例，-alg指定计算算法。命令编写完成后，便可通过subprocess.call()函数将命令执行，分别生成包含研究区域坡度坡向因子的栅格文件。

(3)经度纬度因子。

由于无人机航测的DEM数据中包含研究区的经纬度范围，因此使用GDAL库中的GetGeoTransform()函数即可获取研究区域的经纬度范围，但获取到的坐标是投影坐标，还需要编写函数将投影坐标转换为经纬度坐标，函数主要用到的是TransformPoint()方法。最终得到研究区域的经纬度范围：东经105°46′～105°50′；北纬37°38′～37°43′。基于经度纬度的区间，使用Numpy的linspace()函数分别生成该地区的经度纬度分布矩阵。

4.2 温度推算模型的实现

温度推算模型的实现需要通过GDAL库将各个因子数据的GeoTiff文件进行读取，转化为可以进行运算的矩阵。主要操作如下：

(1)叠加因子权重系数：由于每个因子对该区域的温度影响作用不同，因此需要对每个因子叠加不同权重系数；

(2)叠加总系数：每个因子叠加完权重系数后，推算出的温度区间可能与实际温度分布区间不符，因此需要设置总系数对温度分布区间进行调整；

(3)叠加基准站温度T：前两项运算得出的是目标点与基准站温度的差值，因此需要与基准站温度进行叠加，从而实现对整个区域的温度进行推算。

以上三种对矩阵的操作皆可通过Numpy库进行实现，这里需要注意的是，在执行完每一步操作以后，都要使用Matplotlib库中的imshow方法将运算完成的矩阵以图形化展示，这样可以保证在运算出现异常值时及时发现，同时可以对每一步运算结果进行验证，提升温度推算工作的鲁棒性。

基于用户体验需求，工具在使用时需要隐藏具体的实现细节，同时又需要将权重系数配置接口开放给用户。因此引入了Python的ConfigParser模块，通过此模块可以实现读取配置文件中的权重系数实现温度的推算，因此用户仅需在配置文件中配置好每个因子的权重系数即可生成研究区域的温度分布图。

4.3 温度分布图的生成

基于Python技术对提取到的5个因子通过温度推算模型进行运算以后，最终生成包含每个网格点温度值的GeoTiff文件。虽然基于python的Matplotlib模块可以将温度推算中的每个因子进行图形化展示，方便检查数据的准确性，但是在实际应用中发现Matplotlib模块在绘制最终正果温度分布图时存在很多不足的地方，如无法设置图例样式，无法添加指南针比例尺等。因此决定采用ArcMap制图工具[16-17]制作最终的正果温度分布图。ArcMap是一款强大的制图工具，可以方便地添加图例、比例尺等组件，同时可以灵活地根据数值范围设置不同的配色。

5 实验结果

最终生成的正果温度分布如图5所示，由图中不同标注的区域可以看到，研究区域的西侧区域温度较低，而研究区域北部区域及东南部区域温度偏高，其中温度最高的区域为东北部区域，温度达到了15.3 ℃以上。

图5 正果温度分布

为了对项目建立的温度推算模型进行验证，分别在太阳梁试验区选取三个观测点，对观测点的温度进行逐时监测，统计了2018年3月1日至31日的每小时温度。具体观测位置如图6所示。

图6 监测点位置分布

监测点A位于试验区的北部区域，监测点B位于试验区的西南部区域，监测点C位于试验区的东南部区域。通过统计得出了26日至30日的平均温度：观测点A的平均温度为15.1℃，观测点B的平均温度为15.13℃，观测点C的平均温度为14.9℃。经过与模型推算结果对比，得出推算结果与地区实际温度分布状况比较吻合，达到了预期的目标。同时，温度推算设定的五个因子中，海拔对温度的推算影响最大，其次是坡度和坡向。可见在小区域网格温度推算中，经度和纬度的影响权重较小，基本上可以忽略不计，这一点与大区域的网格点温度推算不同。

6 结束语

基于对温度推算模型及开源栅格处理GDAL库的研究，设计并实现了一种基于DEM数据的小网格温度推算模型，并基于开源栅格处理库GDAL库及科学运算库Numpy库进行实现。首先介绍了研究区概况及模型的设计，其次介绍了实现模型时用到的关键技术，并详细介绍了温度推算模型的实现步骤，最后对温度推算结果进行了分析与验证。基于开源GDAL库实现的温度推算模型相较于传统的实现方式具有可扩展性强且运算速度快等特点。该成果同样可推广应用于其他模型的推算工作，如基于该模型中的坡度、坡向、海拔因子，可以应用于小区域风场分布情况的推算，这也是该项目接下来的研究方向。