GeoTiffCOG icon indicating copy to clipboard operation
GeoTiffCOG copied to clipboard

Bilinear Interpolation

Open anoldfellow opened this issue 1 year ago • 0 comments

Hello,

as I understand, the getElevation code always takes the sample next up-left to the given coordinate. This is no problem for flat areas. But what along a steep coast or in regions like the european Alps? I did some simple experiments with this code. As a result, I calculate the elevation from all four points around the point of interest. The resulting height simply is taken from a bilinear interpolation along x and y path. My code is as follows (not tested well, needs optimyzation):

public float GetElevationAtLatLon(double latitude, double longitude) { float[] h = new float[4]; double h12, h34, xratio, yratio; Utils.GetXYFromLatLon(metadata, latitude, longitude, out int x, out int y);

if (x < 0) throw new Exception($"Error longitude: {longitude} not valid to [{metadata.PhysicalStartLon}, {metadata.PhysicalEndLon}] offsetX: {x}"); if (y < 0) throw new Exception($"Error latitude: {latitude} not valid to [{metadata.PhysicalStartLat}, {metadata.PhysicalEndLat}] offsetY: {y}");

h[0] = GetElevationAtPoint(x, y); h[1] = GetElevationAtPoint(x + 1, y); h[2] = GetElevationAtPoint(x, y - 1); h[3] = GetElevationAtPoint(x + 1, y - 1); // now do bilinear interpolation xratio = (longitude - metadata.PhysicalStartLon - x*metadata.pixelSizeX) / metadata.pixelSizeX; yratio = (latitude - metadata.PhysicalStartLat + (metadata.Height-y-1) * metadata.pixelSizeY) / metadata.pixelSizeY; h12 = (double)h[0] + ((double)(h[1] - h[0])) * xratio; h34 = (double)h[2] + ((double)(h[3] - h[2])) * xratio; return (float) (h12 + (h34 - h12) * yratio); }

anoldfellow avatar Jul 22 '22 13:07 anoldfellow