往期回顾

在前期的数据准备中我们获得了ANUSPLIN插值格式的excel文件,需要进一步通过SPSS转化为固定ASCII码格式的dat文件,并记录下dat文件各表头的格式,准备完成后就可以进行脚本的编制以及后期的数据处理。

因此本文的主要内容分为三个部分:

(1)Excel格式转化为dat格式;

(2)插值脚本的编制及修改说明;

(3)插值数据后期处理。

1. 打开spss后,加载进去excel数据;

2. 设置好对应的表头信息,如下图所示:

需要注意的是,站点的数据类型为字符型,设置好后,另存为固定ASCII码格式的:

在保存的时候注意spss会有个输出的日志文件,下图中标红的部分需要在脚本中用到。

设置好后,可以开始写插值的脚本了,脚本主要用到splina.exe和lapgrd.exe.两个程序。

Splina的脚本文件如下:

0 #插值的单位,0表示无单位

2 #自变量的个数,这里指经纬度

1 #协变量的个数,这里指dem,注意也可以将高程设置为自变量,不用协变量

.1 6206316.3 0 1   ##前两数是对应dem文件中经度的范围(在arcgis中的source看),但要比dem中的范围要大一点才行

.01 5937788.2 0 1 ##前两数是对应dem文件中经度的范围,但也要比dem中的范围要大一点才行

-163.422 7371.14 1 1         ##前两数dem文件的上下限

365 #需要插值的个数,这是是对每天进行插值,故是365个

ssd2007.dat #插值的dat文件名

1000 #站点的个数,可以设置的比实际站点数量大些

5 #标签的宽度,我们的站点都是五位数的,所以是5

(a12,f12.4,f12.4,f12.2,365f12.1)  #数据的格式,与上面spss的输出日志保持一致

ssd2007.res

ssd2007.opt

ssd2007.sur

ssd2007.lis

ssd2007.cov

(注意这里要留一行)

里面的#标注部分根据自己的数据来进行修改即可,```这个符号这里只是用来囊括插值代码的,在写入文件的时候不用写进去。将上述信息写在一个txt文件中,然后另存为一个cmd格式的文件。假设这里另存的文件名称为splina_ssd2007.cmd。Lapgrd的脚本文件如下:

ssd2007.sur

0 #表示输出所有的插值面

.25165 6206316.25165 8000 #前两个是dem文件中的经度范围,8000:插值的分辨率

.1524 5937788.1524 8000 #前两个是dem中的纬度范围,8000:插值的分辨率,最好和dem的分辨率一致

chinadem.txt  #ascii码的dem文件名称

ssd2007_001.grd

ssd2007_002.grd

ssd2007_003.grd

ssd2007_004.grd

ssd2007_005.grd

ssd2007_006.grd

ssd2007_007.grd

ssd2007_008.grd

ssd2007_009.grd

ssd2007_365.grd

#注意这里有个空行才行

同样将上述信息写在一个txt文件中,然后另存为一个cmd格式的文件。假设这里另存的文件名称为lapgrd_ssd2007.cmd。新建一个运行的脚本,命名为run.cmd

splinassd2007s.log

lapgrdssd2007l.log

将代码复制到run.cmd中,双击run.cmd就能够运行得到结果,注意:上述文件必须都放在同一个文件夹下。插值结果是grd格式的,需要在arcgis中批量转换为tif格式。

(1)基于python的grd格式批量转为tif格式

grd格式是传统的arcinfo的格式,不适合进一步的计算,需要将grd转换为tif格式,本文采用arcgis中自带的python来进行批量转换。

具体代码如下所示:

import arcpy

arcpy.env.workspace="D:\\chazhi\\"  #存放数据的文件夹

a=arcpy.ListRasters("*","grd")  #得到文件夹下所有的grd名称

for i in a:

arcpy.RasterToOtherFormat_conversion(i,"D:\\chazhi\\tif\\","TIFF")

将上述代码复制到自带的python上后运行即可。

(2)基于matlab的气象要素异常值处理

插值完后会发现气象要素如降水会存在着负值,负值是由于引进了协变量DEM造成的,可以认为出现的负值的地方降水为0的,尤其是在进入日尺度上插值时,当有n多个这样文件需要处理时,通过arcgis中的栅格计算器来一个个实现尤为费时,本文提供一个基于matlab的处理,同时也能够加入投影信息。

[a,R]=geotiffread('F:\项目\dem.tif');%先导入投影信息

info=geotiffinfo('F:\项目\dem.tif');

e=dir('*.tif');

for i=1:size(e,1)

data=importdata(e(i).name);

data(data<0)=0;

geotiffwrite(e(i).name,data,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

在运行matlab时,需要将工作目录调整到刚刚转换成tif文件的目录下。

以上就是“ANUSPLIN插值专题系列”的全部内容了。接下来,对本专题系列进行一个简要的总结:

(1)气象数据的收集与整理,整理成ANUSPLIN插值所需格式,中间需要用到Excel、SPSS和Matlab,其中Matlab针对长时间序列复杂的气象数据整理;

(2)对气象站点经纬度的转化,将各个站点的经纬度转换成m为单位,需要用到AIbes投影的DEM文件;

(3)DEM文件转化为ASCII码;

(4)脚本文件的构造与修改;

(5)后期插值数据的进一步处理。

欢迎小伙伴们批评指正插值过程中出现的错误,后期会收集各种错误,针对各个错误进行详细的解答。

下期预告:考研资讯|近三年研究生就业率Top10

责任编辑:杨俊凤 尹礼唱 卢雅焱

推文审核:张天舒 梁龙武 骆丹云

总审核:学术无界顾问团

往期回顾: