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

Question

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.

Answer 1

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

Answer 2

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.

See also http://docs.opencv.org/2.4/doc/tutorials/imgproc/imgtrans/warp_affine/warp_affine.html

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

// Machine EPSILON
#define DBL_EPSILON 2.2204460492503131e-16       

// Your transformation
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");
   }

}