针对tif格式DEM数据,由于其数据组成为等间距横向和纵向分布的离散高程点构成,可以看做为普通影像的像素值,点与点之间间距为一个像素,那么对DEM数据的内插,就可以当做在每个像素内均匀增加点,每个点的高程值则由周围四个点通过双线性内插计算结果来确定。
假设当点P(x,y)位于P1(x1,y1)、P2(x2,y2)、P3(x3,y3)和P4(x4,y4)四个点构成的矩阵内,如图所示:
设定四个点高程值为p1、p2、p3、p4,中间P点高程为p,P点高程通过P1-P4四个点高程进行拟合。这里采用距离加权的方法对内插点高程进行赋值,即通过计算P点坐标与其他几个点距离来确定每个值所占权重,假设P点刚好位于矩阵中心,那么p刚好等于p1-p4的均值。
问题的核心在于如何计算P点与周围四个点距离,由于P点在矩阵内,假设一个矩阵边长为1,x和y与周围点坐标的差值都小于等于1,那么可以通过如下公式来计算p点高程值。
通过以上公式,当P点位于矩阵中心时,x-x1=0.5,y-y1=0.5,P1点权重为0.5*0.5=0.25,P2到P4点权重也一样是0.25,如果P点刚好位于P1-P4点任意一点上,通过公式计算结果也是p1-p4。
核心代码:
double result = p1 * ((1 - Math.abs(x - x1)) * (1 - Math.abs(y - y1)))
+ p2 * ((1 - Math.abs(x - x2)) * (1 - Math.abs(y - y2)))
+ p3 * ((1 - Math.abs(x - x3)) * (1 - Math.abs(y - y3)))
+ p4 * ((1 - Math.abs(x - x4)) * (1 - Math.abs(y - y4)));
假设P1 = 100.25, p2= 120.23, p3 = 103.49, p4 = 105.66,P点坐标为(0.2,0.3),计算结果p=112.217,四个点均值为107.4075,P点坐标最接近P2点,P2点最高,因此P点高程大于均值也更符合预期。
具体实现:
public void interpolation(File inputDEM,File outputDEM,int interpolation) {
DemData data = new DemData();
try {
gdal.AllRegister();
Dataset dataset = gdal.Open(inputDEM.toString(), gdalconstConstants.GA_ReadOnly);
Band band = dataset.GetRasterBand(1);
int width = dataset.GetRasterXSize();
int height = dataset.GetRasterYSize();
data.setWidth(width);
data.setHeight(height);
double[] values = new double[width * height];
double[] buf = new double[width];
for (int i = 0; i < height; i++) {
// 一行一行读取
band.ReadRaster(0, i, width, 1, buf);
for (int j = 0; j < width; j++) {
int index = i * width + j;
values[index] = buf[j];
}
}
data.setValues(values);
Driver driver = gdal.GetDriverByName("GTiff");
int destWidth = (width - 1) * interpolation;
int destHeight = (height - 1) * interpolation;
Dataset demDataset = driver.Create(outputDEM.toString(), destWidth, destHeight, 1);
Band demBand = demDataset.GetRasterBand(1);
demBand.SetMetadata(band.GetMetadata_Dict());
double[] resultLine = new double[destWidth];
for (int i = 0; i < destHeight; i++) {
for (int j = 0; j < destWidth; j++) {
double x = (double)j / interpolation;
double y = (double)i / interpolation;
int x1 = (int)x;
int x2 = x1;
int x3 = x1 + 1;
int x4 = x3;
int y1 = (int)y + 1;
int y2 = y1 - 1;
int y3 = y2;
int y4 = y1;
double p1 = data.getValue(y1, x1);
double p2 = data.getValue(y2, x2);
double p3 = data.getValue(y3, x3);
double p4 = data.getValue(y4, x4);
double result = caculateElevation(x, y, x1, x2, x3, x4, y1, y2, y3, y4, p1, p2, p3, p4);
resultLine[j] = result;
}
demBand.WriteRaster(0, i, destWidth, 1, resultLine);
}
demDataset.SetMetadata(dataset.GetMetadata_Dict());
dataset.delete();
demDataset.delete();
} catch (Exception e) {
e.printStackTrace();
}
}
public double caculateElevation(double x, double y, int x1, int x2, int x3, int x4, int y1, int y2, int y3, int y4,
double p1, double p2, double p3, double p4) {
double result = p1 * ((1 - Math.abs(x - x1)) * (1 - Math.abs(y - y1)))
+ p2 * ((1 - Math.abs(x - x2)) * (1 - Math.abs(y - y2)))
+ p3 * ((1 - Math.abs(x - x3)) * (1 - Math.abs(y - y3)))
+ p4 * ((1 - Math.abs(x - x4)) * (1 - Math.abs(y - y4)));
return result;
}