DEM挖填方分析--基于水平参考面计算

基于水平参考高程面计算挖填方比较简单，水平参考高程也就是某一个高程值，只需要计算同一位置上DEM高程到参考面高程即可，挖填方分析实际上就是计算DEM数据与参考面构成的封闭体体积。

DEM数据本身是由一系列等间距横向和纵向分布的高程点构成，其数据组成形式与一般图像类似，相当于每个图像坐标上由像素值变成高程值。一般tif格式DEM数据还附带一个tfw文件用以描述起始点坐标、横向和纵向间隔距离等信息，通过读取tfw文件可以更方便的计算DEM数据。

由于DEM数据为离散的高程数据，为了得到更精确的计算结果，需要对DEM数据进行内插，在DEM内部内插出更多的高程点形成更小的高程格网，以每个格网内高程点作为附近区域平均高程与参考面高程计算差值，并计算体积。当差值为正，此时需要将高出参考面部分挖去，差值即为需要挖掉的深度，通过差值与高程点所在格网面积计算的体积即为该格网处的挖方；当差值为负，此处高程低于参考面，需要进行填方计算，通过同样方式计算填方量。

整个实现过程主要分为三步：

第一步，读取DEM数据，获取每一个栅格的高程值，读取tfw文件，获取起始坐标和每一个栅格长度；

Tfw文件内容：

第二步，对DEM高程值进行内插计算，将每一个栅格分割为n*n的格网，通过双线性内插计算得到格网上每个点的高程值；

第三步，将每个格网上的高程值与参考面高程计算差值，通过单个格网面积计算得到挖填方结果。

具体实现：

public DigFillResultData analysisByElevation(File currentDEM,double referenceElevation,int interpolation) throws Exception{

		// 读取数据
		File currentDEMTfw = new File(currentDEM.toString().replace("tif", "tfw"));
		TfwInfo currentDEMInfo = new TfwReader().read(currentDEMTfw);
		DemData currentData = new DemReader().read(currentDEM);
		// 进行计算
		DigFillResultData result = computeByElevation( currentDEMInfo, currentData,referenceElevation, interpolation);
		return result;
	}
	
private DigFillResultData computeByElevation(TfwInfo currentDEMInfo,DemData currentData,double referenceElevation,int interpolation) {

		double maxDigDeep = 0.0;
		double maxFillDeep = 0.0;
		double digVolume = 0.0;
		double fillVolume = 0.0;
		int width = currentData.getWidth();
		int height = currentData.getHeight();
		DigFillResultData result = new DigFillResultData();
		for (int i = 0; i < height - 1; i++) {
			for (int j = 0; j < width - 1; j++) {
				double currentXCoordinate = currentDEMInfo.getxCoordinate() + j * currentDEMInfo.getxStep();
				double currentYCoordinate = currentDEMInfo.getyCoordinate() + i * currentDEMInfo.getyStep();

				double interXStep = currentDEMInfo.getxStep() / interpolation;
				double interYStep = currentDEMInfo.getyStep() / interpolation;
				for (int x = 0; x < interpolation; x++) {
					for (int y = 0; y < interpolation; y++) {
						double currentXInterCoord = currentXCoordinate + x * interXStep;
						double currentYInterCoord = currentYCoordinate + y * interYStep;

						double currentValue = new ElevationInterpolation().caculateElevation(currentXInterCoord,
								currentYInterCoord, currentData, currentDEMInfo);
						if (currentValue != -32767.0) {
							double difference = currentValue - referenceElevation;
							if (difference > 0) {
								if (difference > maxDigDeep) {
									maxDigDeep = difference;
								}
								digVolume += difference * interXStep * Math.abs(interYStep);
							}
							if (difference < 0) {
								if (Math.abs(difference) > maxFillDeep) {
									maxFillDeep = Math.abs(difference);
								}
								fillVolume += Math.abs(difference) * interXStep * Math.abs(interYStep);
							}
						}
					}
				}
			}
		}
		result.setDigVolume(format(digVolume));
		result.setFillVolume(format(fillVolume));
		result.setMaxDigDeep(format(maxDigDeep));
		result.setMaxFillDeep(format(maxFillDeep));
		return result;
	}