# PostGIS raster - [x,y] from lat/lon

I have coordinates in [lat, lon]. How can I convert them to the rasters "pixel" position [x, y]? I know raster top-left corner, scale and skew.

• Martin, can you tell more on what you want to achieve with this? Calculating the pixel position is not very difficult but a bit ugly and there might be better ways to solve your problem.
– tilt
Jan 30, 2016 at 19:10

Use ST_WorldToRasterCoord(), ST_WorldToRasterCoordX() and/or ST_WorldToRasterCoordY().

In PostGIS you have so called raster editors to set the transformation parameters.

http://postgis.net/docs/manual-2.1/RT_reference.html#Raster_Editors

ST_SetUpperLeft, ST_SetScale, ST_SetSkew...

And raster accessors to transform the numbers.

http://postgis.net/docs/manual-2.1/RT_reference.html#Raster_Accessors

ST_WorldToRasterCoord

The calculation works like this (a small example in C, should work in other languages too). The array notation (order of the params) of the affine transformation params is GDAL conform.

#include <stdio.h>
#include <stdlib.h>

// Machine EPSILON
#define DBL_EPSILON 2.2204460492503131e-16

double trfm[6];

trfm[0] = 7.812837;   // top left x
trfm[1] = 2.3791E-7;  // scaleX is w-e pixel resolution
trfm[2] = 0;          // skew X
trfm[3] = 54.291764   // top left y
trfm[4] = 0           // skew Y
trfm[5] = -2.3791E-7; // scale Y n-s pixel resolution (negative value)

// Transformation from World to Pixel Coodinates
int calcWorldToPixel(double * trfm,
double x, double y,
long * col , long * row) {

double div = (trfm[2]*trfm[4]-trfm[1]*trfm[5]);
if (div<DBL_EPSILON*2) return 0;
double dcol = -(trfm[2]*(trfm[3]-y)+trfm[5]*x-trfm[0]*trfm[5])/div;
double drow =  (trfm[1]*(trfm[3]-y)+trfm[4]*x-trfm[0]*trfm[4])/div;
*col = round(dcol); *row = round(drow);
return 1;
}

int main(void) {
long col =0; long row = 0;
double lon =  7.8128374;
double lat = 54.2917643;
int res = calcWorldToPixel(trfm, lon, lat , &col, &row);
if (res) {
printf(" LON: %f LAT: %f COL: %ld ROW:%ld \n",
lon, lat, col, row);
} else {
printf("Sorry Error\n");
}

}
• This might be a correct solution but the title of the question suggests that it should be done in postgis raster.
– tilt
Jan 30, 2016 at 19:12
• This condition is there why? if (div<DBL_EPSILON*2) return 0; If it is for zero division, this is ot correct, because it will filter out negative div Jan 31, 2016 at 8:32
• I use thr perl style to get an idea if something in the operation is wrong or not. If the operation is valid the func return 1, if not it return 0. By this technique you can formulate things like func1 && func2 && func3 ... Jan 31, 2016 at 17:06