Форумы GIS-Lab.info

Есть желание сделать небольшое приложение, которое на основе данных SRTM будет показывать пользователю высоту в заданной им точке (широта, долгота). Никогда раньше с SRTM не работал. Может быть кто пробовал уже? Какие методы интерполяции лучше использовать (линейную, сплайнами, Делоне или что-то еще) ?

зачем интерполяция?

Ну, как ни крути, не всегда произвольная точка попадет в i,j компоненту матрицы. Поэтому нужно каким-то образом "найти" значение в произвольной точке на основе имеющихся

OMAAGAAD писал(а):Ну, как ни крути, не всегда произвольная точка попадет в i,j компоненту матрицы. Поэтому нужно каким-то образом "найти" значение в произвольной точке на основе имеющихся

Прежде чем пользоваться какой-то моделью, следует выяснить, как она построена. SRTM30 представляет не высоту рельефа в точке, а среднюю высоту по площадке.

Выбор метода интерполирования определяется свойствами функции поверхности, построенной на SRTM. Если она должна быть непрерывна вместе с первой производной, достаточно бикубической интерполяции. Если первая производная не волнует, подойдёт билинейная. Но это всё излишества, поскольку Вам просто нужно значение, соответствующее точности модели. Используйте ближайшее соседство, а лишнюю энергию направьте на что-нибудь полезное.

OMAAGAAD писал(а): не всегда произвольная точка попадет в i,j компоненту матрицы.

всегда, если это не void-ы из первых версий SRTM и они заполнены в производных продуктах типа CGIAR SRTM

http://www.paulinternet.nl/?page=bicubic
C#:

Код: Выделить всё

#pragma strict

using System.Collections;

public class BicubicInterpolator {
	
	//coefficients
	private float
		a00, a01, a02, a03,
		a10, a11, a12, a13,
		a20, a21, a22, a23,
		a30, a31, a32, a33;

	float[,] inData;

	// constructor
	public BicubicInterpolator(float[,] data, int x, int y) {
		float[,] d = new float[4, 4];
		d[0,0] = data[x - 1, y - 1];
		d[0,1] = data[x - 0, y - 1];
		d[0,2] = data[x + 1, y - 1];
		d[0,3] = data[x + 2, y - 1];
		
		d[1,0] = data[x - 1, y - 0];
		d[1,1] = data[x - 0, y - 0];
		d[1,2] = data[x + 1, y - 0];
		d[1,3] = data[x + 2, y - 0];
		
		d[2,0] = data[x - 1, y + 1];
		d[2,1] = data[x - 0, y + 1];
		d[2,2] = data[x + 1, y + 1];
		d[2,3] = data[x + 2, y + 1];
		
		d[3,0] = data[x - 1, y + 2];
		d[3,1] = data[x - 0, y + 2];
		d[3,2] = data[x + 1, y + 2];
		d[3,3] = data[x + 2, y + 2];
		updateCoefficients(d);
	}
	private void updateCoefficients(float[,] d){
		a00 = d[0, 0];
		a01 = -.50f * d[1, 0] +  .50f * d[1, 2];
		a02 = d[1, 0] - 2.50f * d[1, 1] + 2.00f * d[1, 2] -  .50f * d[1, 3];
		a03 =  .50f * d[1, 0] + 1.50f * d[1, 1] - 1.50f * d[1, 2] +  .50f * d[1, 3];

		a10 =  .50f * d[0, 1] +  .50f * d[2, 1];
		a11 =  .25f * d[0, 0] -  .25f * d[0, 2] -  .25f * d[2, 0] +  .25f * d[2, 2];
		a12 = -.50f * d[0, 0] + 1.25f * d[0, 1] - d[0, 2] +  .25f * d[0, 3]
			+  .50f * d[2, 0] - 1.25f * d[2, 1] + d[2, 2] -  .25f * d[2, 3];
		a13 =  .25f * d[0, 0] -  .75f * d[0, 1] +  .75f * d[0, 2] -  .25f * d[0, 3]
			-  .25f * d[2, 0] +  .75f * d[2, 1] -  .75f * d[2, 2] +  .25f * d[2, 3];

		a20 = d[0, 1] - 2.50f * d[1, 1] + 2.00f * d[2, 1] -  .50f * d[3, 1];
		a21 = -.50f * d[0, 0] +  .50f * d[0, 2] + 1.25f * d[1, 0] - 1.25f * d[1, 2]
			- d[2, 0] + d[2, 2] +  .25f * d[3, 0] -  .25f * d[3, 2];
		a22 = d[0, 0] - 2.50f * d[0, 1] + 2.00f * d[0, 2] -  .50f * d[0, 3] - 2.50f * d[1, 0]
			+ 6.25f * d[1, 1] - 5.00f * d[1, 2] + 1.25f * d[1, 3] + 2.00f * d[2, 0]
			- 5.00f * d[2, 1] + 4.00f * d[2, 2] - d[2, 3] -  .50f * d[3, 0]
			+ 1.25f * d[3, 1] - d[3, 2] + 0.25f * d[3, 3];
		a23 = -.50f * d[0, 0] + 1.50f * d[0, 1] - 1.50f * d[0, 2] +  .50f * d[0, 3]
			+ 1.25f * d[1, 0] - 3.75f * d[1, 1] + 3.75f * d[1, 2]
			- 1.25f * d[1, 3] - d[2, 0] + 3.00f * d[2, 1] - 3.00f * d[2, 2] + d[2, 3]
			+  .25f * d[3, 0] -  .75f * d[3, 1] +  .75f * d[3, 2] -  .25f * d[3, 3];

		a30 = -.50f * d[0, 1] + 1.50f * d[1, 1] - 1.50f * d[2, 1] +  .50f * d[3, 1];
		a31 =  .25f * d[0, 0] -  .25f * d[0, 2] -  .75f * d[1, 0] +  .75f * d[1, 2]
			+  .75f * d[2, 0] -  .75f * d[2, 2] -  .25f * d[3, 0] +  .25f * d[3, 2];
		a32 = -.50f * d[0, 0] + 1.25f * d[0, 1] - d[0, 2] +  .25f * d[0, 3]
			+ 1.50f * d[1, 0] - 3.75f * d[1, 1] + 3.00f * d[1, 2] -  .75f * d[1, 3]
			- 1.50f * d[2, 0] + 3.75f * d[2, 1] - 3.00f * d[2, 2] +  .75f * d[2, 3]
			+  .50f * d[3, 0] - 1.25f * d[3, 1] + d[3, 2] -  .25f * d[3, 3];
		a33 =  .25f * d[0, 0] -  .75f * d[0, 1] +  .75f * d[0, 2] -  .25f * d[0, 3]
			-  .75f * d[1, 0] + 2.25f * d[1, 1] - 2.25f * d[1, 2] +  .75f * d[1, 3]
			+  .75f * d[2, 0] - 2.25f * d[2, 1] + 2.25f * d[2, 2] -  .75f * d[2, 3]
			-  .25f * d[3, 0] +  .75f * d[3, 1] -  .75f * d[3, 2] +  .25f * d[3, 3];
	}
	
	public float getValue (int x, int y) {
            //Кешируем возвещение во 2 и 3 степени.
            float x2 = x * x;
            float x3 = x2 * x;
            float y2 = y * y;
            float y3 = y2 * y;

            return (a00 + a01 * y + a02 * y2 + a03 * y3) +
                   (a10 + a11 * y + a12 * y2 + a13 * y3) * x +
                   (a20 + a21 * y + a22 * y2 + a23 * y3) * x2 +
                   (a30 + a31 * y + a32 * y2 + a33 * y3) * x3;
    }
}

Форумы GIS-Lab.info

Интерполяция данных SRTM

Интерполяция данных SRTM

Re: Интерполяция данных SRTM

Re: Интерполяция данных SRTM

Re: Интерполяция данных SRTM

Re: Интерполяция данных SRTM

Re: Интерполяция данных SRTM