3

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 '16 at 19:10
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");
   }

}
  • This might be a correct solution but the title of the question suggests that it should be done in postgis raster. – tilt Jan 30 '16 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 – Martin Perry Jan 31 '16 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 ... – huckfinn Jan 31 '16 at 17:06
2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.