[Gdal-dev] Find elevation data at specific point
Frank Warmerdam
warmerdam at pobox.com
Fri Jun 4 13:46:25 EDT 2004
Gerald Brandt wrote:
> Thanks! What would the pixelx, liney formula be if the image is not
> north up?
>
> I'm working with DEM and DTED and have not yet seen anything that isn't
> north up, but that doesn't mean there aren't any!
Gerald,
Ummm, well it gets more complicated.
I normally invert the geotransform matrix, and then use it in format fashion.
The following code will invert the geotransform[]:
/************************************************************************/
/* InvGeoTransform() */
/* */
/* Invert a standard 3x2 "GeoTransform" style matrix with an */
/* implicit [1 0 0] final row. */
/************************************************************************/
int InvGeoTransform( double *gt_in, double *gt_out )
{
double det, inv_det;
/* we assume a 3rd row that is [1 0 0] */
/* Compute determinate */
det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
if( fabs(det) < 0.000000000000001 )
return 0;
inv_det = 1.0 / det;
/* compute adjoint, and devide by determinate */
gt_out[1] = gt_in[5] * inv_det;
gt_out[4] = -gt_in[4] * inv_det;
gt_out[2] = -gt_in[2] * inv_det;
gt_out[5] = gt_in[1] * inv_det;
gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
return 1;
}
Then you could compute the pixel/line from geo_x, geo_y like this:
pixel = gt[0] + geo_x * gt[1] + geo_y * gt[2]
line = gt[3] + geo_x * gt[4] + geo_y * gt[5]
Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent
More information about the Gdal-dev
mailing list