开源GIS-第六篇:基于openlayer的克里金插值

1.首先获取目标图层点数据

2.基于已有数据训练参数

const variogrm=kriging.train(values,logs,lats,krigingModel,krigingSigma2,krigingAlpha)

3.设置目标区域外包矩形，并在矩形类按照像素大小设置格网

const grid=kriging.grid(polygons,variogram,pixelSize);

4.给每个网格赋值

kinging.plot(canvas,grid,extent,colors);

5,如若有行政边界线对预测结果进行剪裁

kinging.plot(canvas,grid,extent,colors,coordinates);

// Plotting on the DOM
  kriging.plot = function(canvas, grid, xlim, ylim, colors, coordinates) {
	// Clear screen
	var ctx = canvas.getContext("2d");
	ctx.clearRect(0, 0, canvas.width, canvas.height);

	// Starting boundaries
	var range = [xlim[1]-xlim[0], ylim[1]-ylim[0], grid.zlim[1]-grid.zlim[0]];
	var i, j, x, y, z;
	var n = grid.length;
	var m = grid[0].length;
	var wx = Math.ceil(grid.width*canvas.width/(xlim[1]-xlim[0]));
	var wy = Math.ceil(grid.width*canvas.height/(ylim[1]-ylim[0]));
	for(i=0;i<n;i++)
	    for(j=0;j<m;j++) {
		if(grid[i][j]==undefined) continue;
		x = canvas.width*(i*grid.width+grid.xlim[0]-xlim[0])/range[0];
		y = canvas.height*(1-(j*grid.width+grid.ylim[0]-ylim[0])/range[1]);
		z = (grid[i][j]-grid.zlim[0])/range[2];
		if(z<0.0) z = 0.0;
		if(z>1.0) z = 1.0;
		//验证多边形
        var inside = false;
        var xxx=i*grid.width+grid.xlim[0];
        var yyy=j*grid.width+grid.ylim[0];
        // 依次遍历多边形的每个边
        for (var ii = 0, jj = coordinates.length - 1; ii < coordinates.length; jj = ii++) {
          var xi = coordinates[ii][0], yi = coordinates[ii][1];
          var xj = coordinates[jj][0], yj = coordinates[jj][1];

          var intersect =
            ((yi > yyy) != (yj > yyy)) // 判断该点纵坐标是否在线段最高点和最低点之间[注1]
            && (xxx < (xj - xi) * (yyy - yi) / (yj - yi) + xi); // 判断该点向x轴正方向发出的射线是否穿过线段[注2]

          if (intersect) inside = !inside;
        }
    if(inside)
    {
		  ctx.fillStyle = colors[Math.floor((colors.length-1)*z)];
    }
    else
    {
      ctx.fillStyle='#ffffff00';
    }
		ctx.fillRect(Math.round(x-wx/2), Math.round(y-wy/2), wx, wy);
	    }
    };

ps：所调用js参考kinging.js